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SUMMARY 


The  principal  objective  of  this  investigation  was  to  experimentally 
determine  the  feasibility  of  using  acoustic  pulse  echo  techniques  to  mea¬ 
sure  the  thickness  of  sea  ice.  Field  tests  performed  on  both  winter  and 
annual  sea  ice  showed  that  a  pulse  echo  technique  is  entirely  feasible  if 
1)  adequate  coupling  between  the  sound  source  and  the  ice  surface  is  ob¬ 
tained,  and  2)  either  the  source  levels  are  kept  sufficiently  high  or  the 
frequency  is  kept  low  enough  that  the  wavelength  of  the  sound  is  large 
compared  to  the  size  of  the  sound-scattering  centers  in  the  ice.  Our  mea¬ 
surements  of  the  compressional  and  shear  wave  velocities  parallel  to  the 
axis  of  the  ice  indicate  that  such  thickness  measurements  can  be  accurate 
to  within  1.5%. 

The  tests  indicated  that  the  acoustic  impedance  mismatch  at  the  ice- 
water  interface  is  sufficient  to  allow  the  use  of  either  the  compressional 
or  the  shear  wave.  Before  the  shear  wave  could  be  used,  however,  an  effi¬ 
cient  method  would  have  to  be  found  for  establishing  a  reliable  bond  between 
the  ice  and  the  acoustic  source.  We  did  not  obtain  sufficient  data  to 
verify  reports  in  the  literature  that  seasonal  variations  in  both  the  com- 
p^essional  and  shear  wave  velocities  approach  20%. ^ 

Using  the  results  of  this  experiment,  we  empirically  predicted  the 
depth  dependence  of  the  elastic  constants  in  the  ice  and  the  velocity  of 
the  vertically  polarized  shear  wave.  The  average  values  of  these  predic¬ 
tions  are  consistent  with  those  reported  in  the  literature. 

The  attenuation  of  both  the  shear  and  the  compressional  wave  was  then 
calculated  using  a  combination  of  viscoelastic  and  scattering  theory.  The 
theory  was  parametrically  fit  to  measured  values,  and  the  results  compared 
favorably  with  those  of  other  experiments. 

Finally,  preliminary  studies  were  performed  to  determine  to  what 
extent  acoustic  surface  waves  are  propagated  in  sea  ice,  and  to  estimate 
the  feasibility  of  utilizing  these  modes  for  thickness  measurements.  Based 
on  the  results  of  this  study  and  theoretical  considerations,  it  appears 
that  the  symmetric  Lamb  wave  could  be  useful  for  making  ice  thickness 
predictions . 


INTRODUCTION 

Activities  in  ice  covered  regions  require  a  knowledge  of  ice  thick¬ 
ness  to  assure  safety  of  equipment  and  personnel.  This  knowledge  is 
especially  important  during  exercises  that  call  for  the  operation  of 
heavy  equipment,  for  example,  tl.  employment  of  aircraft  on  ice  runways 
for  logistical  support. 

During  the  last  75  years,  many  studies  of  ice  have  been  performed  to 
understand  its  physical,  mechanical  and  chemical  properties.  A  thorough 
discussion  of  this  work  has  been  published  by  Weeks  and  Assur.1  Even 
though  some  of  these  studies  indicate  that  ice  is  a  relatively  good  con¬ 
ductor  of  acoustic  disturbances,  an  efficient  method  of  acoustically  deter¬ 
mining  ice  thickness  from  the  top  surface  still  does  not  exist,  instead, 
indirect  ice  thickness  measurements  are  presently  based  on  other  physical 
characteristics  of  ice  and  water,  such  as  dielectric  permittivity,  dielec¬ 
tric  loss  tangent  in  ice  and  in  seawater,  or  elastic  oscillations  of  the 
total  ice  co\.er. 

Simple  calculations  using  the  idealized  assumption  of  isotropic,  homo¬ 
geneous  ice  indicate  that,  because  of  the  low  acoustic  impedance  mismatch 
at  the  ice-water  interface,  an  acoustic  compressional  wave  incident  on  the 
boundary  will  be  almost  entirely  transmitted  into  the  water  medium.  Con¬ 
versely,  because  the  water  will  not  support  shear  waves,  the  shear  stress 
is  assumed  to  vanish  at  the  boundary  indicating  that  the  wave  will  be 
totally  reflected  at  normal  incidence.  Thus,  neglecting  real  and  apparent 
attenuation,  it  would  appear  from  this  simple  argument  that  the  best  meth¬ 
od  of  determining  ice  thickness  acoustically  from  the  top  surface  would  be 
to  generate  a  shear  wave  in  the  ice  and,  knowing  the  wave  speed,  calculate 
the  thickness  from  the  reflection  of  the  wave  from  the  ice-water  interface. 


This  ideal  model  is,  of  course,  a  gross  oversimplification  of  the 
physics  of  the  problem.  It  is  known  that  sea  ice  is  not  a  simple  medium, 
but  rather  a  highly  complex  inhomogeneous,  anisotropic,  polycrystalline 
viscoelastic  material  for  which  the  acoustic  propagation  characteristics 
may  be  dependent  on  the  growth  and  life  history.  The  principal  areas  of 
concern  to  be  resolved  before  ice  thickness  can  be  reliably  measured  acous¬ 
tically  are  the  velocity  and  attenuation  of  the  various  acoustic  modes  in 
the  medium  and  the  acoustic  impedance  mismatch  of  the  ice-water  interface. 

The  acoustic  velocity  data  reported  for  both  fresh  and  saline  ice 
have  been  acquired  by  one  cf  three  techniques: 

A.  Seismics,2-7  where  the  bulk  sound  velocities  are  determined  indi¬ 
rectly  from  measurements  of  the  longitudinal  plate  waves.  A1 - 
though  this  technique  is  satisfactory  for  "order  of  magnitude" 
estimations,  it  does  contain  several  deficiencies  which  include: 


1.  The  propagation  velocities  of  the  plate  wave  are  highly  dis¬ 
persive  with  the  thickness  of  the  medium  as  was  theoretically 
predicted  by  Lamb,8  Osborn  and  Hart,9  and  Bunney  and  Goodman10 
and  experimentally  verified  by  Bunney  et  al.11  Therefore, 
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unless  the  experimentor  is  very  careful  to  assure  that  the 
wavelength  is  long  compared  to  the  ice  thickness,  large 
variations  in  the  measured  velocity  can  occur. 

2.  The  medium  is  anisotropic,  and  therefore  exhibits  no  true 
"bulk  velocity."  Thus,  measurements  of  long-range  acoustic 
propagation  in  the  ice  layer  yield  information  on  only  one 
component  of  the  velocity  vector,  and  even  this  is  based  on 
the  assumption  that  the  medium  is  homogeneous  which,  of 
course,  it  is  not. 

3.  Because  of  the  distances  involved  with  this  technique,  the 
data  received  at  the  sensor  have  been  both  temporally  and 
spatially  averaged.  The  information  received  therefore  gives 
little  detail  regarding  the  microstructure  of  the  medium. 

B.  Resonance  methods,12*13  where  the  velocities  in  the  medium  are  re¬ 
lated  to  the  normal  modes  of  vibration  of  the  sample.  With  this 
technique  it  must  be  assumed  that  the  test  section  being  measured 
is  homogeneous  because  inconsistencies  in  the  material  will  gen¬ 
erate  dispersive  effects  thereby  altering  the  measured  nodes  and 
antinodes  of  the  resonant  column.  For  natural  sea  ice,  this  is 

an  invalid  assumption. 

C.  Direct  determination  of  the  velocity  by  measuring  the  time  of 
flight  of  an  acoustic  pulse  traversing  ice  samples  of  known  thick¬ 
ness.  Research  using  this  technique  has  been  related  to  velocity 
measurements  in  single  ice  crystals,14'15  measurements  of  lake 
ice  and/or  commercially  available  isotropic  ice,16~1B  or  quasi  in 
situ  velocity  measurements  in  sea  ice.19"21  Of  the  latter,  only 
Bunney  and  Hanse  have  reported  the  internal  velocity  structure  of 
the  medium  and. related  it  to  values  of  previously  measured  bulk 
wave  properties. 

The  attenuation  of  a  sound  wave  propa  rating  in  a  material  that  is  not 
perfectly  elastic,  homogeneous  and  isotroric  can  be  considered  to  be  the 
result  of  two  mechanisms:  (1)  dissipation  processes  originating  from 
internal  friction,  anelastic  behavior  of  the  material,  thermal  dissipation, 
viscous  slippage  at  crystal  boundaries,  etc.,  and  (2)  scattering  originat¬ 
ing  from  the  interaction  of  the  acoustic  wave  with  scattering  centers  in 
the  medium.  The  relative  contribution  of  eacn  of  these  mechanisms  to  the 
total  attenuation  depends  on  the  frequency  of  the  sound  wave.  At  low  fre¬ 
quencies,  the  wavelength  of  the  sound  is  very  large  compared  to  the  scat¬ 
tering  centers;  the  scattering  cross  section,  i.e.,  the  relative  amount  of 
energy  scattered  out  of  the  incident  wave,  is  therefore  extremely  small 
and  the  attenuation  is  due  almost  entirely  to  dissipation.  As  the  fre¬ 
quency  increases,  the  attenuation  due  to  scattering  becomes  more  important 
until,  when  the  wavelength  becomes  approximately  on  the  order  of  the  size 
of  the  scattering  center,  scattering  predominates.  The  attenuation  of 
the  compressional  wave  has  been  experimentally  measured  by  Pounder  and 
Langleben,22  and  the  attenuation  of  both  the  transverse  and  compressional 
waves  in  annual  sea  ice  has  been  mathematically  approximated  by  Hanse  and 
Bunney. 2 3,2 4 
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Although  similarities  exist  between  sea  and  freshwater  ice,  the  dif¬ 
ferences  in  structure  and  composition  are  sufficient  to  preclude  mutual 
substitution  of  experimental  data.  For  a  detailed  comparison  of  fresh¬ 
water  ice  and  sea  ice,  the  reader  is  referred  to  the  excellent  text  on  the 
subject  by  Pounder.19  It  should  also  be  pointed  out  that  Hunkins1*  has 
demonstrated  that  the  physical  properties  of  sea  ice  exhibit  large  seasonal 
variations.  Perennial  ice,  that  which  has  existed  for  more  than  one  sea¬ 
son,  changes  through  the  process  of  free  flooding  during  the  summer  months 
and  differs  greatly  from  annual  ice,  that  which  is  less  than  1  year  old. 


EXPERIMENTAL  TECHNIQUE  AND  RESULTS 

The  experiments  reported  here  were  performed  in  the  Chukchi  Sea  near 
Pt.  Barrow,  Alaska,  during  March  1973.  A  comparative  study  of  the  acoustic 
properties  of  annual  and  perennial  ice  during  this  season  was  planned,  but 
only  annual  ice  samples  were  available  because  of  unavoidable  logistic 
difficulties.  A  comprehensive  investigation  of  the  acoustic  properties  of 
first  year  sea  ice  was  therefore  performed.  The  data  presented  are  from 
that  portion  of  the  study  that  utilized  direct  measurement  techniques  to 
determine  the  depth  dependence  of  the  acoustic  parameters. 


RESULTS  OF  EXPERIMENTAL  MEASUREMENTS 

Numerous  ice  samples  were  acquired  by  extracting  vertical  cores  from 
the  ice  canopy  with  a  10-cm  diameter  SIPRE  coring  tool.  The  average  depth 
of  the  ice  was  approximately  145  cm  and,  with  the  exception  of  one  core 
where  the  skeleton  structure  of  the  ice-water  interface  was  lost,  the 
lengths  of  the  recovered  samples  were  within  2  cm  of  the  measured  ice 
depth.  The  temperature  profile  was  obtained  by  inserting  a  thermometer 
in  small  holes  drilled  along  the  extent  of  the  core  at  6-7  cm  increments. 
The  first  temperature  taken  was  near  the  bottom  of  the  core  because  of  the 
large  temperature  difference  between  that  location  and  the  ambient.  The 
composite*  results  of  the  temperature  tests  for  all  cores  are  given  in 
Figure  1.  It  is  obvious  that  the  technique  employed,  although  used  in 
similar  experiments,16  * 17  is  not  adequate  to  predict,  the  temperature  pro¬ 
file.  Even  after  2  or  3  minutes  the  bottom  of  the  core  cools  by  4-5°. 
Because  15-20  minutes  were  required  to  obtain  an  entire  contour  by  this 
method,  errors  of  10°C  or  more  are  likely.  As  will  be  explained  later, 
the  best  estimate  of  the  core's  temperature  profile  during  the  experiment 
is  given  by 

|  T  =  (0.137Z  -  21.7) UC  ,  (1) 

where  z  is  the  depth  from  the  surface  in  centimeters  and  the  average  ‘ 
ambient  air  temperature  was  -21.7°C.  This  gradient  is  shown  as  the 
broken  line  on  Figure  1. 

•The  data  paints  presented  in  Figures  1  and  3-6  were  obtained  by  plotting 
the  raw  data  as  a  function  of  depth  from  the  surface,  interpolating 
between  points  and  then  averaging  over  all  samples  at  5-cm  increments. 
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Figure  1.  Measured  temperature  profile  vs  depth  from  surface  of  ice 
(lines  extending  from  data  points  represent  data  spread ). 
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After  the  temperature  measurements,  the  cores  were  cut  into  lengths  of 
6  or  7  cm  and  each  segment  was  referenced  relative  to  its  position  in  the 
sample.  The  compressional  wave  velocities  parallel  and  perpendicular  to 
the  sample  axis  were  then  obtained  for  each  section  using  the  instrumen¬ 
tation  diagramed  in  Figure  2.  In  this  system,  a  pulse  timing  generator 
furnishes  pulse  length  and  rate  information  to  a  pulsed  power  oscillator 
and  a  trigger  signal  to  the  oscilloscope.  The  oscillator  in  turn  provides 
a  cw  pulse  of  the  desired  frequency,  rate  and  length  to  a  BaTi03  transducer 
coupled  to  the  sample.  The  transmitted  acoustic  signal  is  received  by  a 
second  transducer  coupled  to  the  opposite  side  of  the  test  section,  ampli¬ 
fied,  and  transmitted  to  the  oscilloscope.  The  time  of  flight  of  the 
acoustic  pulse  in  the  sample  is  determined  by  measuring  with  a  calibrated 
time  delay  the  oscilloscope  sweep  time  from  the  origin  to  the  first  re¬ 
ceived  signal,  and  then  correcting  for  electronic  time  losses  (7  ysec  for 
this  system) .  Knowing  the  specimen  length,  the  velocities  can  then  be 
calculated.  This  method  is  accurate  to  approximately  1-1/2%,  where  the 
largest  error  is  related  to  the  length  measurement. 


TRANSDUCER  TRANSDUCER 


Figure  2.  Black  diagram  of  instrumentation  for  acoustic  velocity 
measurements . 

Figures  3  and  4  show  the  results  of  the  velocity  measurements.  The 
velocity  profiles  exhibited  in  the  graphs  increase  rapidly  for  the  first 
50-55  cm,  and  then  decrease  to  minima  at  ~70  cm  depth.  It  is  unlikely 
that  these  minima  are  characteristic  of  all  annual  sea  ice  because  the 
ice  had  obviously  undergone  severe  "rafting"  earlier  in  the  growth  season. 
Visual  inspection  of  the  cores  revealed  a  "frazil"  layer  in  the  region  of 
70  cm  indicating  that  the  area  had  been  "rafted"  at  least  once.  Below 
the  70-cm  minima  the  profiles  increase  again  to  maxima  at  roughly  100  cm 
and  then  decrease  toward  the  seawater  sound  velocity  at  the  ice-water 
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interface.  The  average  velocities  measured  for  the  vertical  and  horizontal 
directions  were  3840  m/sec  and  3635  m/sec,  respectively.  The  velocities 
compare  favorably  with  results  obtained  by  the  bulk  velocity  measurements 
reported  by  other  investigators  as  shown  in  Table  I,  page  9. 


Figure  3.  Longitudinal  velocity  parallel  to  sample  axis  ae  a  function  of 
depth  from  the  ice  surface  ( lines  extending  from  data  points 
represent  data  spread ). 


7 


Table  I.  Comparison  of  reported  data  for  acoustic  velocities  and  elastic  moduli. 
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Figure  6 .  Salinity  as  a  function  of  depth  from  the  ice  surface  ( lines 
extending  from  data  points  represent  data  spread). 
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The  experimental  results  of  the  sample  density  measurements  are  shown 
in  Figure  5  and  demonstrate  that  for  depths  greater  than  approximately  100 
cm  the  density  increases  markedly,  indicating  a  higher  liquid  saturation. 
From  the  "idealized"  temperature  gradient  shown  in  Figure  1,  the  in  situ 
medium  is  predicted  to  have  a  temperature  of  -8.0°C  at  depths  greater  than 
100  cm.  Pounder19  and  Assur20show  that  as  the  temperature  of  ideal  sea¬ 
water  is  lowered  through  the  freezing  point  a  two-part  system  consisting 
of  solid  H2O  and  uniform  brine  is  formed .  Continuing  to  lower  the  temper¬ 
ature  results  in  the  solidification  of  sodium  sulfate  at  about  -8°C  form¬ 
ing  a  three-part  system  of  solid  H2O,  solid  Na2S04*10Ii20  and  brine.  This 
situation  exists  to  approximately  -22°C  where  sodium  chloride  (NaCl*2H20) 
solidifies.  Thus,  the  increase  in  density  below  depths  of  100  cm  is  due 
to  the  decreasing  amount  of  solid  salt  present. 

After  the  temperature  and  density  measurements,  the  segments  were 
melted  and  the  salinity  was  measured.  The  average  salinity  profile  is 
shown  in  Figure  6.  Although  not  as  large  as  expected,  particularly  near 
the  bottom,  the  depth  contour  is  characteristic,  showing  increased  salini¬ 
ty  near  the  top,  which  is  caused  by  a  thermal  gradient  across  the  elongated 
brine  cells,  and  near  the  bottom  because  of  gravity. 


APPROXIMATION  OF  THE  ELASTIC  MODULI  AND  TRANSVERSE  VELOCITY  PROFILES 

It  was  demonstrated  earlier  that  the  acoustic  parameters  measured 
exhibit  dependence  with  depth  into  the  medium.  Utilizing  the  results  of 
the  reported  experiment,  an  approximation  to  the  structure  and  a  predic¬ 
tion  of  the  average  values  of  the  transverse  velocity  and  elastic  moduli 
can  be  made.  To  accomplish  this  it  is  assumed  that  the  measurement  fre¬ 
quency  is  sufficient  to  ignore  the  viscous  terms  in  the  propagation  equa¬ 
tion,  or  equivalently,  that  the  viscous  terms  do  not  gx'eatly  affect  the 
velocity.  An  elastic  medium  is  thus  assumed,  for  which,  according  to 
elasticity  theory,  the  longitudinal  and  transverse  wave  velocities  are 
given  by 


where 

Cl  is  the  velocity  of  the  longitudinal  wave 
Gp  is  the  velocity  of  the  transverse  wave 

E  is  Young's  modulus 

y  is  the  shear  modulus 

A  is  the  Lame"  constant 
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p  is  the  density  of  the  medium 

a  is  Poisson's  ratio. 


There  are  two  empirical  equations  in  the  literature  that  associate 
Young's  modulus  (E)  with  the  chemical  properties  of  the  medium.  Pounder 
and  Stalinski17  related  this  parameter  to  the  salinity  (S)  by 

E  -  (9.75  -  0.242  S)  x  1010  dyn/cm2  ;  (3) 

however,  their  work  was  performed  at  constant  temperature.  Langleben18 
related  Young's  modulus  to  the  brine  content  (o)  by 

E  »  (10.0  -  35.1  V)  x  1010  dyn/cm2  .  (4) 

The  brine  content  in  this  equation  can  be  determined  from  the  salinity  and 
temperature  by  the  relations  25 

V  «  S  (£?#■  -  2.28  j  .  2.06°  <  0  <  -0.5° 

v  a  s  jr||—  +  0.930^  -  8.2°  <  0  <  -2.06°  (5) 

V  =■  5  ^‘jj—  +  1  •  189^  -22.9°  <  0  <  -8.2°  , 

where 

S  is  the  salinity  in  parts  per  thousand 
v  is  the  brine  content  in  parts  per  thousand 
0  is  the  temperature  in  degrees  Celsius. 

Because  the  temperature  dependence  is  included  in  Eq.  5  it  is  expected  to 
give  a  more  accurate  estimate  of  Young's  modulus  than  Eq.  4. 

As  was  discussed  earlier,  the  method  used  to  measure  the  temperature 
profile  during  the  experiment  could  lead  to  large  errors  because  of  the 
temperature  difference  between  the  sample  and  ambient,  and  the  time  re¬ 
quired  to  perform  the  tests.  To  reconstruct  the  in  situ  conditions,  we 
have  referred  to  Untersteiner 's26  approximation  of  the  thermal  conductiv¬ 
ity  in  sea  ice, 

Ki  Ki,f  T-273  *  C6r 

where 

S(z)  ■  salinity  at  depth  z 
0  ■  0,28  cal  cm2/g  sec 
“  0.00486  cal/cm  sec  °K 
T  ■  temperature  in  degrees  Kelvin. 
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Examination  shows  that  fox  absolute  temperature  values  greater  than  a  few 
degrees  and  salinities  on  the  order  of  a  few  parts  per  thousand,  the  sec¬ 
ond  term  of  this  equation  can  be  neglected,  yielding  a  constant  thermal 
conductivity.  These  conditions  are  applicable  to  the  problem  being  dis¬ 
cussed,  and  allow  the  assumption  of  a  linear  temperature  gradient  between 
the  average  ambient  during  the  experiment  (-21.7°C)  and  the  seawater 
freezing  temperature  (-1.8°C). 

Thus,  using  the  relationship 

T  »  (0.137z  -  21.7) °C 

and  the  experimental  values  of  the  longitudinal  velocity,  density  and 
salinity,  the  depth  dependence  of  the  brine  content  can  be  approximated 
along  with  Young's  modulus,  Poisson's  ratio,  the  transverse  velocity,  the 
shear  modulus  and  the  Lame"'  constant.  The  results  of  these  calculations 
are  shown  in  Figures  7  through  12.* 


*Each  of  the  results  presented  can  be  reduced  to  an  empirical  relationship 
with  depth.  However,  the  temptation  to  do  so  has  been  resisted  because 
the  data  may  be  unique  to  the  medium  for  the  1973  season.  The  dependency 
with  depth  shown  in  the  figures  has  been  derived  from  the  raw  data 
acquired  and  thus  will  show  structure  which  may  not  be  reproducible  in 
future  experiments. 
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Figure  9.  Predicted  values  of  Poisson's  ratio  vs  depth  from  surface . 
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Figure  11,  Predicted  values  of  the  shear  modulus  vs  depth  from  surface. 


DISCUSSION  OF  THE  RESULTS  OF  THE  PROFILE  MEASUREMENTS 


In  order  to  evaluate  the  credibility  of  the  order  of  magnitude  esti¬ 
mates,  the  results  of  our  predictions  and  those  reported  by  others  are 
compared  in  Table  I.  Wide  variations  do  occur  for  some  of  the  values,  but 
in  general  the  predicted  values  and  the  average  of  the  experimental  data 
compare  very  well  with  those  reported  in  the  literature. 

Shear  wave  velocity  profile  experiments  have  yet  to  be  successfully 
performed  in  sea  ice  because  of  the  difficulty  in  establishing  a  reliable 
bond  between  the  ice  and  shear  wave  generating  transducers.  This  diffi¬ 
culty  is  largely  because  the  brine  "leaches"  and  produces  a  highly  saline, 
liquid  surface  layer.  Washing  the  bonding  surface  with  fresh  water  has 
been  partially  successful  in  laboratory  and  in  situ  tests.  For  example, 
the  experiment  shown  in  Figure  13  was  attempted  during  the  field  tests 
reported  here.  Channels  in  the  ice  canopy  were  cut  to  a  depth  of  -45  cm 
and  located  a  known  distance  apart.  The  sides  of  the  channel  were  washed 
with  fresh  water  to  eliminate  the  saline  boundary  layer.  The  shear-sensi¬ 
tive  transducers  were  quickly  coupled  to  the  ice  on  diametrically  opposite 
sides  of  the  test  section  at  a  depth  of  15  cm.  The  electronics  used  were 
the  same  as  shown  in  Figure  2.  The  time  of  flight  and  distance  were  mea¬ 
sured,  after  which  the  sample  size  was  incrementally  decreased  and  the 
process  repeated. 

The  average  value  of  this  measurement  was  1923  m/sec  and  is  shown  as 
a  triangle  in  Figure  10.  Further  laboratory  measurements  of  shear  wave 
velocity  on  sea  ice  samples  show  an  average  of  1870  m/sec,  compared  to  a 
predicted  average  profile  velocity  of  1874  m/sec.  Because  sea  ice  is 
anisotropic,  it  should  be  pointed  out  that  these  values  are  for  the  verti¬ 
cally  polarized  shear  wave.  The  horizontally  polarized  wave  was  measured 
in  a  separate  experiment  and  will  be  discussed  later  in  this  report. 

Figure  7  reinforces  the  earlier  discussion  on  density  by  showing 
that  the  brine  content  increases  rapidly  for  depths  greater  than  100  cm. 
However,  unlike  the  "ideal"  seawater  case  where  the  sodium  sulfate  solid¬ 
ifies  discontinuous ly  at  approximately  -8.0°C  this  calculation  shows  that 
at  least  some  of  the  salts  solidify  uniformly  with  decreasing  temperature. 


TRANSMITTER  RECEIVER 


SEA  ICE 
SURFACE 


Figure  13.  Cross  section  of  vn  situ  shear  wave  experiment. 


Although  Young's  modulus  (Figure  8)  decreases  precipitously  below 
-8°C,  indicating  a  rapid  decrease  in  the  material  elasticity,  the  form  of 
the  shear  modulus  (Figure  11)  is  physically  more  interesting.  Because 
this  parameter  must  vanish  at  the  ice-water  interface,  and  because  in  the 
calculation  it  appears  to  be  uniformly  approaching  that  asymptote,  it  is 
not  anticipated  that  a  significant  acoustic  impedance  mismatch  for  the 
shear  wave  will  occur  at  the  boundary.  If  the  calculation  were  definitive, 
the  interface  would  not  be  a  good  acoustic  reflector  for  the  transverse 
mode.  However,  since  strong  echoes  of  this  mode  have  been  observed  in 
situ,  it  is  obvious  that  greater  consideration  must  be  given  to  the  thoory 
for  the  transition  region. 

PREDICTIONS  OF  ACOUSTIC  ATTENUATION  IN  SEA  ICE 

Before  the  feasibility  of  measuring  sea  ice  thickness  acoustically 
can  be  determined,  the  attenuation  of  sound  in  the  medium  as  a  function 
of  frequency  must  be  either  measured  or  estimated.  As  discussed  earlier, 
the  attenuation  of  a  sound  wave  propagating  in  a  material  that  is  not  per¬ 
fectly  elastic,  homogeneous  and  isotropic  is  largely  dependent  on  two 
mechanisms:  (1)  dissipation  processes  originating  from  internal  friction, 

anelastic  behavior  of  the  material,  thermal  dissipation,  viscous  slippage 
at  crystal  boundaries,  etc.,  and  (2)  scattering  originating  from  the 
interaction  of  the  acoustic  wave  with  scattering  centers  in  the  medium. 


LONG  WAVELENGTH  ISOTROPIC  SOLUTION 

For  wavelengths  that  are  long  compared  to  the  size  of  the  scattering 
center,  the  calculation  of  acoustic  velocity  and  attenuation  can  be  de¬ 
rived  ignoring  the  contributions  due  to  scattering.  Even  though  not  ail 
of  the  internal  mechanisms  contributing  to  the  dissipation  processes  are 
known,  the  general  theory  can  be  developed  by  grouping  all  attenuation 
into  a  common  source.  This  is  accomplished  mathematically27  by  replacing 
the  shear  modulus  (M)  and  the  Lams'*  constant  (A)  in  the  elastic  stress - 
strain  relationship  with  the  first  order  differential  operators 

M  "  V  +  V'  & 

5  (’ 

A  -  A  ♦  A'*  £  . 

where  the  unprimed  terms  denote  the  elastic  and  the  primed  terms  denote 
the  viscous  (or  attenuative)  contributions. 


Consider  an  acoustic  shear  wave*  moving  in  the  positive  [downward)  z- 
direction.  The  displacement  potential  is  given  by 


where 


iji  *  Ae1  -  az) 


(8) 


u)  is  the  circular  frequency 
t  is  the  time 

a  is  the  complex  wavenumber. 

This  displacement  potential  obeys  the  wave  equation 


at2 


M  V2  i|>  , 


(9) 


where  p  is  the  density  of  the  medium  and  V2  is  the  Laplacian  operator. 
Substituting  Eqs.  7  and  8  into  Eq.  9  and  carrying  out  the  operations  gives 
the  relation 


pur 


a2  (p  +  iwiO 


CIO) 


where  u'  is  the  dissipative  coefficient.  For  an  attenuative  medium,  the 
complex  wavenumber  can  be  written  as 


a  *  k  -  it 


Cil) 


where 

k  is  the  running  vector,  u/c 
c  is  the  wave  velocity 
T  is  the  generalized  attenuation. 

Making  this  substitution  into  Eq.  11  and  separating  the  real  and 
imaginary  components  yield  the  simultaneous  equations 

pw2  =»  y(k2  -  t2)  +  2u"ujkT 

0  ■  uip'Ck2  -  t2)  -  2pkr  . 


(12) 


Setting 


r-H_ 


*The  derivation  for  the  compressional  wave  propagation  is  similar  to  that 
presented  here,  the  difference  being  the  introduction  of  the  (A  +  2M) 
operator  in  the  wave  equation,  Eq.  9,  rather  than  the  shear  modulus  (M) 
operator.  The  displacement  potential  must,  of  course,  be  associated  with 
the  compressional  rather  than  the  transverse  wave. 


these  equations  have  the  solutions 


*  ■  LW1  ♦  rV  -1! 
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y2  ,  1  POJ2  (Vi  +  P.V  ♦  1  | 

2  V  \  l+  ^2  j 

where  both  y  and  R  may  he  frequency  dependent.  The  procedure  from  this 
point  is  to  introduce  a  mathematical  model  of  the  medium  composed  of 
elastic  and  viscous  constants  that  are  frequency  independent,  and  then 
fit  the  results  to  experimental  data. 

Maxwell*8  suggested  that  viscoelastic  materials  could  be  represented 
by  an.  elastic  element  (spring)  in  series  with  a  viscous  element  (dashpot) 
as  shown  in  Figure  14.  Voigt29  placed  these  elements  in  parallel  as  in 
Figure  15.  Detailed  mathematical  analyses  of  these  as  well  as  more  com¬ 
plex  models  are  covered  in  standard  texts30-32  and  will  not  be  discussed 
here.  However,  it  can  be  ^hcwn  that  for  the  Maxwell  model 
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Figure  14.  Maxwell's  mechanical  representation 
of  visaoelaetic  solids. 
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Figure  15.  Voigt's  mechanical  representa¬ 
tion  of  viscoelastic  solids. 


For  the  Voigt  model 


where 


V  "  Ey 

Rw  «  (Ev/nv)cj  . 


(17) 


(18) 


Very  few  solids  behave  like  either  the  Maxwell  or  the  Voigt  model. 
However,  since  more  complicated  models  become  extremely  involved  mathe¬ 
matically,  and  because  models  specifying  a  greater  number  of  parameters 


require  more  experimental  evidence  to  substantiate  their  validity,  using 
a  single  Maxwell  or  Voigt  element  is  a  convenient  method  of  obtaining  a 
first  order  approximation  of  a  viscoelastic  solid's  mechanical  properties. 
This  is  especially  true  if  there  are  only  limited  data  on  the  material, 
as  is  the  case  for  sea  ice. 


SHORT  WAVELENGTH  SOLUTION 

When  sea  ice  freezes,  the  brine  in  the  water  concentrates  in  small 
pockets.  The  result  is  frozen  water  with  imbedded  scattering  centers 
filled  with  a  highly  saline  brine.  Although  these  inclusions  can  have 
various  geometries,  the  author  has  observed  that  the  vast  majority  in 
annual  sea  ice  are  spheroidal.  This  may  not  hold  for  other  types  of  sea 
ice,  e.g.,  multi-year  ice,  but  the  acoustic  and  microscopic  analysis  nec¬ 
essary  for  evaluation  of  other  types  of  sea  ice  has  not  been  done.  For 
the  purposes  of  this  calculation,  it  will  be  assumed  that  all  of  the 
scattering  centers  are  brine-filled  spheres. 

To  evaluate  the  attenuation  due  to  scattering,  distribution  functions 
in  the  medium  must  be  either  known  or  assumed.  Because  the  necessary  data 
are  presently  unavailable,  and  to  minimize  the  complexity  of  the  calculation, 
the  following  assumptions  have  been  made: 

(a)  No  interaction  exists  between  the  scattering  centers,  and  the 
wave,  once  scattered,  will  not  be  re-scattered  back  into  the 
sound  beam.  Since  second  order  scattering  should  have  only  a 
small  effect,  particularly  in  the  Rayleigh  scattering  region, 
this  assumption  should  not  greatly  affect  the  calculation 
accuracy. 

(b)  There  are  N  scattering  centers  per  unit  volume  and  the  centers 
are  uniformly  distributed  throughout  the  medium.  From  obser¬ 
vations  in  the  laboratory,  this  assumption  appears  to  be 
reasonably  valid. 

(c)  All  of  the  scattering  centers  are  the  same  size,  with  a  radius 
a0.  In  fact,  the  brine  pocket  sizes  vary  widely;  however,  if 
radius  a0  is  selected  as  the  average  bubble  size,  then  it  can 
be  argued  that  only  the  larger  bubbles  in  the  distribution  will 
significantly  affect  the  scattering. 

Using  these  assumptions ,  the  attenuation  due  to  scattering  can  be 
written 

^scattering  "  <T  ^ao^  * 


where  Y(a0)  is  the  scattering  cross  section  of  a  single  scatterer  of 
radius  a0  (see  Appendix  A  for  the  cross  section  calculation) .  Because 
of  (c),  above,  this  equation  predicts  the  minimum  attenuation  to  be 
expected  from  tiis  mechanism. 
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Laboratory  measurements  of  scattering  centers  in  annual  sea  ice  give 
radii  on  the  order  of  0.065  cm  and  concentrations  on  the  order  of  25  holes 
per  cubic  centimeter.  Using  the  results  of  the  cross  section  calculation 
in  Appendix  A  and  the  relation 


ka  * 


(20) 


where  v  is  the  frequency  in  hertz,  c  is  the  velocity  in  centimeters/second, 
and  a  is  the  radius  in  centimeters,  the  attenuation  due  to  scattering  for 
values  of  ka<l  and  ka>15  can  be  approximated  by 


a 
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and 


COmp  2.294  x  10'3  +  3.253  x  10"3  V4 
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where  v  is  in  megahertz.  For  values  of  l<ka<15,  the  approximations  of 
Eqs.  21  and  22  are  still  more  or  less  valid  for  order  of  magnitude  cal¬ 
culations,  but  will  invariably  yield  predictions  that  are  significantly 
low. 


TOTAL  ATTENUATION  APPROXIMATION 

As  was  discussed  earlier,  the  attenuation  in  sea  ice  is  composed  of 
two  parts, 

a  m  dissipation  +  “scattering  •  (2^J 

Before  the  total  attenuation  can  be  estimated,  the  attenuation  due  to 
dissipation  must  be  considered.  The  experiment  shown  in  Figure  13  gives 
the  results  necessary  for  an  approximation  to  the  solution  of  the  two- 
element  models.  During  this  experiment,  the  total  attenuation  and  ve¬ 
locity  were  measured  in  situ  15  cm  below  the  surface  in  annual  sea  ice. 

At  100  kHz,  the  results  of  these  measurements  gave  a  longitudinal  wave 
velocity  of  3790  m/sec,  a  vertically  polarized  shear  wave  velocity  of 
1923  m/sec,  an  average  density  of  0.937  g/cm1,  and  a  generalized  atten¬ 
uation  of  0.13  dB/cm.  Using  these  values,  and  assuming  the  elastic  and 
viscous  constants  to  be  frequency  independent,  both  the  Voigt  and  the 
Maxwell  models  were  solved  as  a  function  of  frequency.  The  results  are 
shown  in  Figure  16,  along  with  the  scattering  attenuation  for  the  shear 
wave.  From  this  representation,  it  can  be  concluded  that  for  frequencies 
above  approximately  5  kHz  the  Maxwell  model  fails  because  in  nature  there 
is  no  observed  frequency  region  over  which  the  attenuation  remains  con¬ 
stant  .  Thus,  assuming  the  Voigt  solution,  the  total  attenuation  can  be 
approximated  by 
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and  V  is  in  megahertz. 


The  solution  to  Eq.  24  is  shown  in  Figure  17  along  with  some  exper¬ 
imental  data  of  Langleben  and  Pounder22  on  compressional  wave  attenuation 
taken  in  situ  in  Tanquary  Fiord,  Ellesmere  Island,  during  1967.  Although 
it  appears  that  the  theory  predicts  values  of  attenuation  several  decibels 
too  high  in  the  .1-5  kHz  region,  the  form  of  the  frequency  dependence  is 
consistent  in  that  the  form  of  the  observed  relaxation  is  predicted. 
Further,  the  fact  that  the  theory  is  quite  sensitive  to  the  radius  of  the 
scattering  centers  could  account  for  this  difference  between  the  prediction 
and  the  observed  data.  Differences  in  brine  pocket  sizes  and  hence  devi¬ 
ations  in  the  attenuation  could  be  expected  in  any  event,  since  Langleben 
measured  the  attenuation  on  2 -year  sea  ice  in  a  fiord,  while  the  curves 
were  empirically  fit  to  data  obtained  on  annual  ice  in  the  Arctic  Ocean. 


DEPTH  DEPENDENCE  OF  THE  ATTENUATION 


The  recent  experimental  measurements  by  the  author21  have  demonstrated 
that  for  annual  sea  ice  the  velocity  of  longitudinal  waves  exhibits  de¬ 
pendence  on  the  depth  into  the  medium  because  of  variations  in  density, 
salinity,  and  temperature.  From  these  data  a  shear  wave  velocity  profile 
was  calculated  which  predicted  a  continuous  decrease  with  depth. 


Unfortunately,  sufficient  data  do  not  exist  to  similarly  predict  the 
depth  dependence  for  attenuation.  For  this  determination,  the  form  of 
both  the  elastic  and  viscous  coefficients  must  be  known.  However,  if  a 
near-elastic  medium  is  assumed,  it  can  be  shown  that  for  a  Voigt  solid 


pnar 

2C* 


(26) 


29 


"It 

% 


Therefore,  provided  the  depth  dependence  of  the  viscosity  (n)  does  not 
decrease  faster  than  C3,  the  attenuation  will  increase  with  depth. 


ACOUSTIC  REFLECTIVITY  AT  THE  ICE -WATER  INTERFACE 

In  considering  the  feasibility  of  using  acoustic  techniques  to  measure 
ice  thickness,  a  major  factor  will  be  the  amount  of  energy  reflected  from 
the  ice-water  interface.  For  an  acoustic  wave  incident  on  a  boundary  be¬ 
tween  two  ideal  media,  the  relative  energy  distribution  calculation  is 
straightforward  and  is  presented  in  detail  in  Appendix  A.  The  resulting 
equations  were  solved  using  average  values  for  the  density  and  the  com- 
pressional  and  transverse  velocity  of  sea  ice  as  determined  in  a  previous 
section  of  this  report.  The  solution  for  the  compressional  and  transverse 
waves  is  shown  in  Figures  18  and  19,  respectively. 

Neglecting  attenuation  effects,  this  simple  model  indicates  that  at 
normal  or  near  normal  incidence  the  most  efficient  method  of  acoustically 
determining  ice  thickness  would  be  through  the  implementation  of  the 
transverse  mode  because  of  unit  reflectivity  at  the  interface.  There  are 
at  least  two  difficulties  with  this  approach: 

(a)  The  true  velocity  profile  of  the  transverse  wave  is  unknown. 

It  has,  however,  been  predicted  to  be  highly  dependent  upon 
the  amount  of  liquid  brine  present  in  the  medium — the  amount 
varies  as  a  function  of  temperature.  Further,  a  very  sharp 
transition  of  the  medium  occurs  in  the  lower  few  centimeters 
of  the  ice,  i.e.,  the  skeleton  structure.  Therefore,  without 
a  great  deal  more  information,  which  is  presently  unavailable, 
it  is  impossible  to  ascertain  the  acoustic  impedance  mismatch 
at  the  interface  and  thus  predict  the  reflected  energy  in  any 
given  acoustic  mode. 

(b)  An  efficient  method  of  introducing  the  transverse  wave  into 
the  medium  has  yet  to  be  devised.  The  standard  technique  of 
coupling  shear  wave  generating  transducers  onto  the  surface 
fails  for  the  case  of  sea  ice,  principally  because  the  saline 
layer  generated  by  the  brine  leaching  processes  precludes  a 
good  bond  between  the  transducer  and  the  medium. 

Several  methods  were  attempted  to  introduce  sufficient  acoustic  en¬ 
ergy  into  annual  sea  ice  to  observe  reflections  from  the  interface.  The 
only  successful  technique,  wherein  an  acoustic  pulse  is  generated  using 
a  low-level  explosive  charge,  is  shown  in  Figure  20.  A  slit  is  cut  in 
the  ice  a  known  distance  from  the  explosive  source  to  eliminate  the  direct 
path  to  the  sensor,  a  10-kHz  omnidirectional  hydrophone  implanted  in  the 
side  of  the  slit.  The  explosion  excites  all  of  the  acoustic  modes,  and 
these,  in  general,  propagate  at  different  velocities.  The  signal  received 
at  the  sensor  should  therefore  contain  contributions  from  each  mode  that 
are  separated  temporally  according  to  their  different  velocities.  Figures 
21  and  22  show  typical  returns  from  this  experiment.  Because  the  direct 
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Figure  18.  Relative  energy  distribution,  congressional  wave  incident  on 
ice-water  boundary  ( C ^  =  2840  m/sea ;  p  =  0.94  g/cc ). 


Figure  19.  Relative  energy  distribution,  transverse  wave  incident  on  ice- 
'water  boundary  (CT  =  3840  m/secs  C-,  =  1870  m/eec ;  P  *  0.94 
g/cc).  L  r 
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path  has  been  interrupted,  the  first  arrival  at  the  sensor  is  the  compres- 
sional  wave  echo  from  the  interface  and  is  indicated  by  arrow  1  in  the 
figures.  The  initial  signal  is  followed  by  reverberation  until  the 
arrival  of  the  shear  wave  reflection  denoted  by  arrow  2.* 


Figure  20.  Diagram  of  reflectivity  experiment. 


Correcting  for  delays  in  the  oscilloscope  triggering  circuit  and 
averaging  over  all  measurements,  the  shear  wave  velocity  was  1746  m/sec. 
This  value  is  125  m/sec  lower  than  that  predicted  by  the  velocity  calcu¬ 
lations.  However,  there  are  two  shear  modes  that  can  be  excited,  ver¬ 
tical  and  horizontal,  each  of  which  can  propagate  at  different  speeds 
because  of  the  anisotropic  nature  of  the  medium.  It  is  impossible  to 
determine  which  mode  has  been  observed  in  the  experiment,  but  it  is  be¬ 
lieved  that,  because  of  the  nature  of  the  sound  source,  the  horizontally 
polarized  mode  is  the  principal  contributor.  The  results  of  the  in  situ 
tests  are  therefore  not  necessarily  inconsistent  with  the  predicted 
values . 

Figure  23  shows  the  result  of  a  similar  experiment  performed  on 
summer  sea  ice.  Even  though  this  medium  is  quite  different  from  annual 
sea  ice,  the  form  of  the  acoustic  echo  from  the  interface  is  very  simi¬ 
lar.  The  compressional  and  shear  wave  contributions  are  designated  by 
arrows  1  and  respectively. 

There  are  two  obvious  differences  in  the  structure  of  the  return 
echoes  for  summer  and  annual  sea  ice  that  can  be  correlated  to  differ¬ 
ences  in  the  medium.  First,  the  acoustic  returns  in  summer  ice  do  not 
exhibit  the  reverberation  observed  for  the  annual  ice.  This  indicates 
that  the  medium  is  more  homogeneous  during  the  summer  months,  which,  in 
fact,  it  is  because  of  the  free  flooded  conditions.  Second,  the  well" 
defined  return  echoes,  particularly  for  the  shear  wave,  are  indicative 
of  an  abrupt  transition  at  the  interface;  i.e,,  the  skeleton  layer  is 
not  well  formed  during  the  summer  months. 


^Because  of  triggering  difficulties,  T=0  in  Figures  21-23  does  not 
represent  the  initiation  of  the  sound  source. 


Sweep  =0.5  msec/di v,  Amp  =0.01  V/div 
Ice  depth  =  1 .5  m 

Source  to  receiver  distance  =  3.6  m 

Figure  21.  Typical  return  from  pulse  echo  experiment  in  annual  sea  ice. 


1  2 


Sweep  =0.5  msec/div,  Amp  =  0.005  V/div 

Ice  depth  =  1.5  m 

Source  to  receiver  distance  =  5  m 

Figure  22.  Typical  return  from  pulse  echo  experiment  in  annual  sea  ioe. 


1  2 

Sweep  =  1  msec/div,  Amp  =  1  V/div 

Ice  depth  =  2,286  m 

Source  to  receiver  distance  =  2  m 

Figure  23.  Typical  return  from  pulse  echo  experiment  in  summer  sea  ice. 

Unfortunately,  certain  difficulties  during  the  summer  ice  experiments 
did  not  allow  a  calculation  of  the  trigger  delay  so  that  the  average 
velocities  cannot  be  directly  determined.  However,  assuming  the  compres- 
sional  velocity  is  3840  m/sec,  the  transverse  wave  velocity  is  calculated 
to  be  1765  m/sec.  Therefore,  although  it  has  been  shown4  that  the  acous¬ 
tic  parameters  exhibit  some  seasonal  dependence,  this  effect  has  not 
necessarily  been  observed  during  the  experiments  reported  here. 

Although  the  experiments  performed  have  not  directly  measured  the 
acoustic  impedance  mismatch,  they  have  demonstrated  that  sufficient  energy 
is  reflected  from  the  interface  to  monitor  the  thickness  using  acoustic 
pulse  techniques,  at  least  for  annual  and  summer  sea  ice. 

SEISMIC  MEASUREMENTS 

As  discussed  in  the  introduction,  many  seismic  studies  have  been  per¬ 
formed  to  measure  acoustic  velocities  in  sea  ice.  Great  significance  has 
been  placed  on  the  speed  of  the  plate  wave  to  determine  the  corapressional 
wave  velocity,  and  the  frequency  of  the  air-coupled  wave4*6  »33  to  determine 
ice  thickness.  However,  other  acoustic  modes  such  as  the  Lamb  waves8  and 
the  Rayleigh  wave34  could  be  as  important  for  the  determination  of  these 
parameters  as  the  plate  or  air-coupled  wave.  Some  preliminary  experimen¬ 
tation  was  performed  to  better  understand  the  mechanism  of  long-range 
acoustic  transmission  parallel  to  the  ice  cover. 

High  frequency  propagation,  i.e.,  wavelengths  much  smaller  than  the 
ice  thickness,  will  not  contribute  to  the  energy  observed  at  significant 


distances  from  the  sound  source.  Thera  are  two  reasons  for  this:  (1) 

The  attenuation  of  the  sound  wave  increases  at  least  as  the  fourth  power 
of  the  inverse  wavelength  so  that  the  energy  is  dissipated  very  quickly 
for  the  frequencies  being  discussed,  and  (2)  the  radical  velocity  profile 
through  the  medium  causes  extensive  refraction  of  the  wave  front.  Using 
the  horizontal  velocity  profile  discussed  in  a  previous  section  of  this 
report  and  the  computer  program  RAY0Q5  given  in  Appendix  B,  the  refraction 
was  calculated  for  a  sound  source  located  on  the  surface,  near  mid -depth 
and  at  the  interface.  The  result  of  this  calculation  is  shown  in  Figures 
24-26.  It  is  obvious  from  this  representation  that  the  only  condition 
under  which  appreciable  energy  will  be  transferred  is  for  the  source  to 
be  located  at  mid-depth  because  of  the  sound  channel  at  that  location. 

Even  this  condition  may  be  unique  because  of  the  formation  history  of 
the  1973  ice  cover. 

To  examine  the  longer  wavelength  contributions  to  the  transmitted 
acoustic  energy,  the  experiment  illustrated  in  Figure  27  was  performed. 

A  line  array  of  geophones  composed  of  eight  elements  (four  sensing  verti¬ 
cal  displacements  and  four  sensing  horizontal  displacements)  was  located 
oil  the  ice  such  that  each  sensing  station  was  a  known  distance  from  the 
source  of  the  acoustic  disturbance  (a  small  explosive  charge).  The  re¬ 
ceived  signal  at  each  of  the  geophones  was  recorded  on  an  analytic 
recorder  for  later  laboratory  analysis.  Typical  results  of  the  horizon¬ 
tally  oriented  array  are  displayed  in  Figure  28.  Results  for  the  verti¬ 
cally  oriented  array  are  shown  in  Figure  29. 

Although  the  ice  surface  was  fairly  rough,  the  geophones  were  lo¬ 
cated  as  equidistantly  as  possible.  Thus,  in  the  representation  used 
in  Figures  28  and  29  the  velocity  of  the  wavefronts  is  given  by  the 
slope  of  the  line  connecting  the  arrival  times  of  initial  impulse  on 
the  individual  traces.  For  the  horizontal  geophone  array,  there  were 
two  obvious  wavefront  arrivals,  one  travelling  at  3050  m/sec  and  one  at 
1638  m/sec.  Except  for  the  air  wave  at  the  17-meter  geophone  there  is 
no  clear  separation  in  the  arrival  times  of  other  acoustic  modes.  This 
could  be  caused  by  the  long  oscillation  times  and  the  close  proximity  of 
the  geophones,  because  on  several  shots  there  were  signal  interferences 
that  indicated  the  possible  arrival  of  another  acoustic  mode.  An  example 
is  shown  by  the  broken  line  on  Figure  28.  A  best  guess  for  this  data 
indicates  that  this  mode  has  a  velocity  of  1445  m/sec,  which  is  very 
close  to  the  velocity  of  sound  in  water.  Future  experimentation  with  the 
sensors  located  further  from  the  sound  source  should  separate  this  wave 
front. 

The  intermediate  velocity  wave  has  bear,  previously  reported  by 
Oliver,3  who  was  unable  to  explain  its  origin.  Given  that  Cg  is  the 
velocity  of  the  Ra>  "Leigh  wave  and  Cp  and  Cg  are  the  velocities  of  the 
transverse  and  longitudinal  wave,  respectively,  Viktorov35  shows  that 
the  Rayleigh  wave  velocity  can  be  determined  from  the  equation 

n3  -  8i)a  +  8(3  -  25*)ti  -  16(1  -  52)  a  0  ,  (27) 


where 
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Figure  24.  Sound  ray  plot  for  annual  eea  ice ;  source  at  surface. 


’ ay  plot  for  annual  sea  ice;  source  at  mid-depth. 


Using  the  known  values  of  the  longitudinal  velocity,  density,  and  plate 
velocity  for  ice,  and  the  identities 


_ 1 _ 

pCl-cr2) 

E  (l-o) 
p  (1+a) (l- 2a) 

E 

2p(l+a) 


(28) 


where  p  is  the  density,  E  is  Young's  modulus  and  a  is  Poisson’s  ratio,  a 
value  of  Op  =*  1763  m/sec  was  derived.  This  is  in  excellent  agreement  with 
the  reflectivity  measurements  discussed  earlier.  Equation  27  was  then 
solved  for  the  Rayleigh  velocity  for  this  medium  and  found  to  be  1677 
m/soc,  a  value  very  near  that  observed  for  the  intermediate  speed  wave. 


Figure  27.  Diagram  of  aeiomio  experiment. 

Other  evidence  indicates  that  the  observed  disturbance  is  due  to  the 
Rayleigh  surface  wave.  It  is  known  that  this  wave  is  propagated  with  an 
elliptical  particle  motion  that  is  unique  to  this  particular  mode.  Thus, 
the  orientation  of  displacement  sensors,  such  as  geophones,  is  unimportant 
in  the  detection  of  this  wave.  Examination  of  Figures  28  and  29  shows 
that  this  disturbance  is  detected  by  both-  the  vertically  and  horizontally 
oriented  geophone  arrays,  while  the  plate  wave  excites  only  the  horizon¬ 
tally  sensitive  geophones.  This  indicates  that  the  particle  motion  for 
this  mode  is  indeed  elliptical.  Finally,  Bunney  and  Goodman36  have  shown 
that  the  Rayleigh  wave,  once  excited,  is  the  predominant  mode  of  energy 
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Result  of  pilot  study  to  measure  Lamb  wares  in  cornual  sea  ice;  vertically  polarized  geophones 
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transfer  and  propagates  with  very  little  attenuation,  which  is  consistent 
with  the  behavior  of  the  observed  wave. 

To  fully  interpret  the  results  of  the  seismic  experiment,  the  theory 
of  Bunney  and  Goodman10  has  been  extended  to  allow  excitations  in  a  plate 
bounded  by  air  on  one  side  and  water  on  the  other  (see  Appendix  A) .  The 
predicted  velocity  dispersion  for  the  case  of  sea  ice,  using  the  acoustic 
parameters  previously  discussed  in  this  report,  is  shown  in  Figure  30. 
There  are  two  surface  wave  modes  predicted:  (1)  the  symmetric  Lamb  wave, 
which  for  small  thickness-to-wavelength  ratios  (k  h<l)  is  indistinguish¬ 
able  from  the  plate  wave  (Cp)  and  appears  to  asymptotically  approach  the 
shear  wave  velocity  (CqO  as  the  layer  thickness  increases,  and  (2]  .he 
Rayleigh  wave,  or  the  antisymmetric  Lamb  wave,  which  appears  to  propagate 
near  the  Rayleigh  wave  velocity.  These  two  waves  were  the  principal 
excitations  observed  during  the  experiment. 


CONCLUSIONS  AND  RECOMMENDATIONS 

The  ability  to  efficiently  and  accurately  determine  the  thickness 
of  sea  ice  using  acoustic  pulse  echo  techniques  is  dependent  on  several 
factors.  The  more  critical  problems,  those  associated  with  the  velocity 
and  attenuation  of  the  acoustic  wave  in  the  medium  and  with  the  acoustic 
impedance  mismatch  at  the  ice-water  interface,  have  been  investigated  in 
the  laboratory  and  in  situ  for  both  summer  and  annual  sea  ice.  The  re¬ 
sults  of  these  measurements  indicate  that,  if  sufficient  coupling  between 
the  sound  source  and  the  medium  can  be  obtained,  and  if  the  acoustic  fre¬ 
quency  is  maintained  such  that  the  wavelength  is  long  compared  to  the 
scattering  centers  in  the  medium,  tho  determination  of  ice  thickness  using 
pulse  echo  methods  is  entirely  feasible.  This  technique  was  demonstrated 
in  situ . 

For  thicknesses  up  to  2.5  meters  (the  maximum  depth  available  during 
the  field  tests]  reflections  of  both  the  compressional  and  transverse 
waves  were  observed.  Although  the  compressional  wave  velocity  measure¬ 
ments  exhibited  radical  fine  structure  with  depth  into  the  medium,  the 
average  values  varied  by  l"ss  than  1.5%.  This  is  indicative  of  the  ac¬ 
curacy  with  which  the  ice  thickness  can  be  measured.  Either  the  compres¬ 
sional  or  the  transverse  wave  can  be  used  for  the  measurements.  Because 
of  the  largo  acoustic  impedance  mismatch  for  the  shear  wave  at  the  ice- 
water  boundary,  the  energy  returned  is  several  decibels  greater  than  that 
for  the  compressional  wave.  This  advantage  may  be  negated,  however,  by 
the  relative  difficulty  in  generating  this  acoustic  mode. 

The  small  explosive  charges  used  as  sound  sources  during  these  experi¬ 
ments  worked  very  well  for  generating  all  of  the  acoustic  modes,  including 
the  sheax*  wave,  Because  of  the  time  required  to  set  up  the  apparatus, 
this  technique  would  be  very  inefficient  for  surveying  large  areas;  thus 
other  methods  of  coupling  energy  into  sea  ice  should  be  studied. 
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Hunkins4  has  reported  seasonal  variations  in  both  the  compressional 
and  shear  wave  velocities  approaching  20%.  We  did  not  obtain  sufficient 
reliable  data  during  these  tests  to  verify  this  result.  Further,  the 
theory  shows  that,  for  high  frequencies,  scattering  processes  are  the 
largest  contributor  to  the  attenuation  of  the  sound  wave  in  the  medium. 

For  the  annual  ice  observed,  the  attenuation  becomes  significant  for  fre¬ 
quencies  above  approximately  50  kHz  and  is  dependent  on  the  size  of  the 
scattering  centers,  in  this  case  the  brine  cells.  Because  the  size  of 
these  cells  can  have  a  large  seasonal  variation,  a  dramatic  effect  may  be 
observed  in  the  attenuation  of  sound  in  the  medium.  It  is  recommended 
that  tests  be  conducted  to  measure  the  seasonal  dependence  of  the  acoustic 
parameters  in  the  medium. 

The  empirical  relationships  developed  from  the  experimental  data  for 
predicting  the  depth  dependence  of  the  acoustic  parameters  (such  as  the 
transverse  wave  velocity)  appear  to  become  invalid  near  the  transition 
layer  between  the  ice  and  seawater.  Because  a  reliable  prediction  of  the 
relative  amount  of  energy  each  mode  reflects  from  this  interface  is  de¬ 
pendent  on  a  thorough  understanding  of  the  acoustic  properties  of  this 
transition  layer,  it  is  recommended  that  a  program  be  initiated  to  study 
the  physics  of  this  "skeleton  layer." 

Preliminary  experiments  have  shown  that  energy  can  be  propagated 
parallel  to  the  surface  of  the  ice  canopy  in  both  the  Rayleigh  and  Lamb 
modes.  Theoretical  results  have  shown  that  the  Lamb  waves  are  highly  dis¬ 
persive  with  the  thickness  of  the  layer.  The  introduction  of  pure  tone 
acoustic  waves  into  the  media  may  well  provide  first  order  approximations 
to  the  averago  ice  thickness  over  significant  ranges.  It  is  recommended 
that  this  technique  receive  further  attention. 
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APPENDIX  A  -  THEORIES 


THEORY  OF  SOUND  WAVE  PROPAGATION  VELOCITY  AND  ATTENUATION  IN  VISCO¬ 
ELASTIC  MEDIA 


A.  Long  Wavelength  Solution 

1.  Derivation  of  the  General  Equations  for  Wavenumber  and 
Attenuation  in  an  Isotropia  Medium 


For  wavelengths  that  are  long  compared  to  the  size  of 
the  scattering  centers,  scattering  theory  can  be  ignored  in 
deriving  the  formulation  of  acoustic  velocity  and  attenuation. 
The  introduction  of  viscoelasticity  is  accomplished  by  the 
replacement  of  the  Lame"'  constants  A  and  M  in  the  elastic 
stress-strain  relations  by  the  first  order  differential 
operators 


A.x.r£ 


(Al) 


where  the  unprimed  terms  are  the  elastic  and  the  primed  terms 
the  viscous  contributions . 

Consider  an  acoustic  shear  wave*  moving  in  the  positive 
(downward)  x-direction.  The  displacement  potential  is  given 
by 

*  .  tJ- (Mt  -  (A2) 

where  w  is  the  circular  frequency,  t  is  the  time,  and  a  is  the 
complex  wavenumber.  The  potential  obeys  the  wave  equation 


p  s  M  V2  *  , 


(A3) 


3t 


where  p  is  the  density  of  the  medium.  Substitution  of  Eqs.  A1 
and  A2  yields 


p  O.  =  !£ 

at2 


•The  derivation  for  compressional  wave  propagation  is  similar  to 
that  presented  here  for  shear  waves,  the  difference  being  the 
introduction  of  A+2M  in  the  wave  equation  rather  than  the  shear 
modulus  M.  The  displacement  potential  must,  of  course,  be 
associated  with  the  compressional  rather  than  the  shear  wave. 
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or 


p  w2  =  a2(p  +  iwp") . 
For  a  dissipative  medium. 


(A5) 


a  =  k  -  ix,  (A6) 

where  k  is  the  running  vector  O/c)  and  T  is  the  attenuation. 
Making  this  substitution  and  separating  real  and  imaginary 
components  yield  the  simultaneous  equations 

pte2  =  u(k2  -  x2)  +  2p "(nkt  rA71 

0  =  wy'Ck2  -  t2)  -  2ykt  .  L  J 
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Setting 


these  equations  have  the  solution 


t  -  -1} 

(A8) 

1.2  _  1  pai2  f-s/l  +  R2u>2  +  1  | 

J  ' 


where  both  p  and  R  may  be  frequency  dependent. 
2.  Mathematical  Models  of  Sea  Ice 


Historically,  there  have  been  four  different  rheological 
models  used  to  represent  sea  ice.  These  models  relate  the 
stress  tensor  as  a  function  of  density,  strain  and  time  rate 
of  strain.  Only  two  models  will  be  discussed  here,  that  of 
Maxwell  and  that  of  Voigt.  The  linear  standard  model  and  the 
four-element  model  are  combinations  of  these  two,  and  their 
introduction  would  only  lead  to  undue  complexity. 

a.  The  Maxwell  Model2* 


Consider  the  mechanical  representation  shown  in 
Figure  Al,  where  the  ice  is  represented  by  a  spring  denot¬ 
ing  the  elastic  element  (1%,)  in  series  with  a  dashpot 
denoting  the  viscous  element  (r^) .  Both  Em  and  nm  are 
frequency  independent.  Maxwell  defined  the  stress-strain 
relationship  for  this  viscoelastic  solid  as 
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(i  t  3) 


(A9) 


where  a  is  the  stress,  S  is  the  strain  and  <S  is  the 
relaxation  time  of  the  material.  If  the  wave  equation 
(Eq.  A3)  is  to  be  applicable,  Hooke's  law  states  that 
the  relationship  between  stress  and  strain  must  be  of 
the  form 


cr  »  2  MS,  (A10) 

which  may  be  accomplished  by  writing  M  as  an  operator. 


(AH) 


or,  in  terms  of  the  elastic  and  viscous  constants  of 
Figure  Al, 


M 


Figure  Al,  Maxwell's  meohaniaal 
representation  of 
visooelastia  solids. 


Thus,  substituting  into  Eq.  A10  and  simplifying  give 


d  m 
-nr  a  +  —  o 

dt  % 


(A12) 
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(A16) 
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b.  The  Voigt  Model2 s 


Figure  A2  shows  the  mechanical  representation  of  the 
Voigt  Model.  In  this  configuration,  a  spring  representing 
the  elasticity  (Ev)  is  in  parallel  with  a  dashpot  repre¬ 
senting  the  viscous  term  (riv)  •  As  in  the  previous  model, 
both  Ev  and  nv  are  assumed  to  be  frequency  independent. 

The  stress -strain  relationship  is  given  by 


0  * 
v 


<V\5t 


(A17) 


Figure  A2.  Voigt's  mechanical 
representation  of 
viaooelastio  solids. 


Using  techniques  and  arguments  similar  to  those  for  the 
Maxwell  model  above,  it  can  be  shown  that  the  solutions  for 
the  Voigt  model  are  given  by  the  equations 
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B.  Extension  of  the  Theory  to  Include  Anisotropy 

Thus  far  in  the  theory,  only  the  solutions  for  an  isotropic 
medium  have  been  discussed.  It  is  known,  however,  that  for  many 
materials,  including  naturally  occurring  ice,  the  acoustic  properties 
exhibit  a  dependence  on  the  direction  of  the  propagation  of  the  wave 
relative  to  the  growth  axis,  i.e.,  the  media  p re  anisotropic.  A 
method  of  including  anisotropy  in  the  theory  follows. 

Stresses  and  strains  in  the  medium  can  be  specified  in  terms 
of  the  resolved  displacements  in  a  unit  cube  of  the  material.38 
There  are  six  components  of  stress  and  six  of  strain,  which  can  be 
written  in  terms  of  tensor  notation  as 

°ij  “d  sij  • 

where  i,j  denotes  the  x-,  y-,  z-  or  1-,  2-,  3-component  of  the  stress 
or  strain.  The  force  on  a  unit  cube  in  the  i^  direction  is  given  by 

3a.  . 

Fi  -  ^ydz  (i J  -  1.2.3)  ,  (A19) 

j 

where  repeated  indices  indicate  a  summation.  The  tensor  strains  are 
defined  by 


(A20) 


where  the  u's  are  the  displacements  of  the  body  along  the  coordinate 
direction  xj . 

Utilizing  the  reduced  tensor  convention  for  the  stresses 
al  ‘  ail;  *2  "  °22;  "  °33 


23 


°32;  a5 


13 


°31’  a6  "  al2  ’  °21 


and,  similarly,  for  the  strains 


S1  '  SIV  S2  -  ®22*  S3  ■  s33 


7  S4  ■  S23  ■  S32>  7  SS  ’  S13  •  S32>  7  S6  *  SI2  -  S21  • 


the  generalized  Hooke's  law  can  then  be  written 
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*i  “  Cij  Sj  (i,j  “  1*2,... 6)  (A21) 

and  it  can  be  shown  that  for  the  elastic  constants 

cij  ■  cji*  CA22> 

There  are  a  maximum  of  21  elastic  constants  for  the  most  unsymraetri- 
cal  crystal.  As  the  degree  of  symmetry  increases,  the  number  of 
elastic  constants  decreases  until  there  are  only  two  for  an  isotropic 
medium.  These  are  the  Lame'*  constants  X  and  u»  where 

x*2v-cn  ■■Cjj.Cjj. 

*  *  ci2  "  C13  "  C23  "  C21  "  C31  =  C32  ’  ^A23'1 

U  "  C44  “  C55  “  C66  ' 

and  all  other  constants  are  zero. 

To  include  viscosity  in  this  theory,  each  of  the  elastic  con¬ 
stants  is  replaced  by  a  differential  operator  of  the  proper  form  for 
the  model  being  considered.  Therefore,  a  complete  solution  will  re¬ 
quire  the  determination  of  two  or  more  elastic  and  viscous  constants 
for  each  operator  used. 

C.  Short  Wavelength  Solution 

For  wavelengths  on  the  order  of  or  shorter  than  the  size  of  the 
scattering  centers  in  the  medium,  the  energy  loss  due  to  scattering 
must  be  considered. 

Assume  a  plane  wave  propagating  through  an  isotropic  solid  to 
be  scattered  by  a  sphere  of  radius  "a."  Both  the  wave  and  the  scat- 
terer  are  characterized  by  their  longitudinal  (ki  and  k2)  and  shear 
(<1  and  «2)  wavenumbers .  The  energy  scattered  out  of  the  sound  beam 
by  the  sphere  is  equal  to  the  time  rate  at  which  energy  is  being  car¬ 
ried  away  across  a  spherical  surface  concentric  with  the  scatterer 
and  of  radius  "b"  >  "a."  The  scattered  energy  is  then  written 
mathematically  as 


1  1 

[  3u  3u  3u 

^scattered  “If 

|  axr  ST  +  °yr  St  *  °zr  3t~ 

A  1 

scattered 


dA  ,  (A24) 
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where  axr,  Cyr,  and  ozr  are  the  complex  stress  components  acting  in 
the  three  rectangular  axes  on  a  surface  normal  to  the  radius  vector 
r  and  where  ux,  Uy,  and  uz  are  ..the  complex  displacements  in  the  same 
three  axes.  Since  the  final  expression  for  the  scattered  energy 
oust  be  real,  and  both  the  fl’s  and  u's  are  generally  complex,  only 
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the  real  part  of  Eq.  A24  will  be  taken.  Rewriting  Eq.  A24  in 
spherical  coordinates  and  determining  the  real  part  gives 


’'scattered 


iw  F  |  la  u*  ♦  a6  u*  +  a  u*l 
~  J  |  [  rr  r  6r  8  zr  z] 


'  [  °rr  ur  +  a0r  U9  +  °zr  uz]  I  2rb2  sin8d6  . 
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The  displacements  of  the  sound  waves  can,  in  general,  be  written 


*  grad  cf>  +  curl  tjl. 


(A26) 


where  is  the  scalar  potential  and  ^  is  the  vector  potential  (0,  0, 
ijj) .  The  displacements  for  the  incident,  scattered  and  internal*  waves 
can  be  written  in  terms  of  Eq.  A26  with  the  incident  wave  potentials 
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the  scattered  wave  potentials 
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and  the  internal  wave  potentials 
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*the  waves  inside  the  scattering  sphere 
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where  jm  hm  and  p^  are  the  spherical  Bessel  function,  the  spherical 
Hankel  function  and  the  Legendre  polynomials,47  respectively.  The 
stress  and  displacement  components  in  terms  of  the  potentials  are39 

C14kV)  (A30) 

a,  -  o 

<j)r 


1  3<b  1  3Z  ,  .. 

r"36  *  r  3r5F 


(A31) 
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Applying  Eqs.  A28,  A30,  and  A31  to  Eq.  A2S  yields  the  scattered 
energy: 


^scattered  -  ^  | b J •  ] 


(A33) 


The  scattering  cross  section  is  defined  as  the  ratio  of  the 
scattered  energy  to  the  total  energy,  or 


mvwwwWa Ww 

Etotal 


(A34) 


For  an  incident  plane  compress ional  wave  with  amplitude  coefficients 

im  -  ir  c-i)B+1  C2o+i) 
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the  total  energy  in  the  absence  of  a  scatterer  is 


PjW3 

Etotal  " 


(A35) 


and  the  scattering  cross  section  is 
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The  final  solution  depends  on  the  amplitudes  Am  and.  Bm  which  must 
be  evaluated  in  terms  of  the  boundary  conditions  defining  the  stress 
and  displacement  at  the  scattering  sphere  surface. 

The  scattering  cross  section  for  a  shear  wave  incident  on  a 
spherical  scatterer  is 


comp 


4fr 


m=o 


2m+l  k 


|Aj2  +  m(m+l)  |  Bffl|  2  j 


(A37) 


where  again  the  boundary  conditions  are  solved  for  the  A^s  and  B^s. 
The  solution  of  this  equation  using  the  best  known  values  for  the 
parameters  shows  that,  except  in  the  region  where  the  ratio  of  the 
circumference  to  the  wavelength  is  approximately  1:1,  the  scattering 
cross  section  obeys  the  relationship 


y  -  - — — ~ "  7  ■  dP/cm2 

Shear  4.151  x  10  ?  +  2.35  x  10  51  V^a'2,2 


v**a 6 


(A38) 


0009  2.163  x  10'7  +  0.017  vV 


where  y  is  the  scattering  cross  section,  a  is  the  scattering  radius, 
and  v  is  the  frequency  in  megahertz. 

To  determine  the  effect  of  scattering  on  the  acoustic  attenua¬ 
tion  of  sea  ice,  the  distribution  of  the  scatterers  in  size  and  space 
must  be  considered. 
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II.  THEORY  OF  ACOUSTIC  REFLECTION  AND  TRANSMISSION  AT  THE  ICE-WATER 

INTERFACE 

The  equations  to  be  developed  here  are  for  isotropic,  homogeneous 
media  whoso  extent  is  very  much  larger  than  the  acoustic  wavelength. 
Complexities  such  as  the  anisotropic  nature  and  the  viscoelastic  proper¬ 
ties  of  sea  ice  are  not  considered,  since  much  of  this  theory  is  already 
complete.40"44 

Elementary  tests  on  acoustics,  such  as  Lindsay,45  show  that  for  an 
acoustic  wave  striking  normal  to  an  interface  such  as  that  shown  in 
Figure  A3  (where  the p  's  are  the  densities,  the  A's  are  amplitudes,  and 
the  C's  are  the  velocities  of  the  acoustic  wave  in  the  media)  the  relative 
amplitude  of  the  reflected  wave  is  given  by 


where  2  =  (o 2C2) / (p jC^)  =  the  acoustic  impedance  of  the  boundary.  The 
above  model  shows,  in  a  very  simplified  way,  that  the  acoustic  energy 
reflected  is  wholly  dependent  upon  the  impedance  mismatch  of  the  interface. 


Figure  A3.  Aoouatio  wave  reflection  and  tranmiasion  at  an  interface  for 
a  normal  incident  wave. 

For  the  case  of  an  incident  wave  not  normal  to  the  boundary,  consider 
the  diagram  of  Figure  -A4,  where  medium  1  (the  ice)  is  assumed  to  be  iso¬ 
tropic  and  homogeneous  and  medium  II  (the  water)  is  assumed  to  be  an  ideal 
fluid.  Then  the  potential  functions  are  given  by 
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Figure  A4.  Acoustict  wave  reflection  and  transmission  at  an  interface  for 
non-normal  incidence . 
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Finally,  the  amplitude  of  the  transmitted  congressional  wave  is 
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From  these  equations,  the  energy  distribution  relative  to  that  of  the 
impinging  wave  can  then  be  calculated  as  a  function  of  incident  angle. 
The  results  of  these  calculations,  using  the  best  known  values  of  the 
necessary  parameters,  are  shown  in  Figures  18  and  19  of  the  text. 


III.  ENERGY  OF  THE  ACOUSTICALLY  EXCITED  SURFACE  WAVES  ON  A  PLATE 

Consider  a  plane  wave  travelling  in  an  ideal  fluid  striking  a  plate 
of  thickness  (h)  at  an  angle  G0.  The  following  mathematical  development 
calculates  the  energy  per  unit  length  per  second  flowing  in  the  y-direction 
(see  Figure  A5),  where  it  is  assumed  that  the  plate  is  bounded  on  each 
side  by  a  (but  not  necessarily  the  same)  fluid.  The  problem  is  to  deter¬ 
mine  if  enhancement  occurs  near  the  Rayleigh  and  Lamb  wave  excitation 
angles . 


Figure  A5.  Schematic  representation  of  the  problem . 
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If  the  displacement  potentials  are  <f>  and  then  the  displacement 


u  ■  grad  <t>  +  curl  $  , 


(A45) 


where  the  potential  of  an  incoming  wave  with  angular  frequency  u),  wave 
number  k0  and  unit  amplitude  is 


i(uit-k  sin  0  y-k  cos  9  x) 
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the  potential  of  the  reflected  wave  is 


i(ut  -  k  sin  0  y  +  k  cos  0  x)  .  . 
=  e  o  o  o  o  ' 


the  potential  cf  the  transmitted  wave  is 


-kFsineF  y  -  kp  cos  0p  x)  , 


and  the  potentials  <j>  and  \p  in  the  elastic  medium  are 
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The  intensity  of  the  wave  in  the  i^  direction  is  given  by37 

I.'  -  ilx(u*  S  +u.  S*  )+  2y(u*  S.  +u  s!  )  >  ,  (A52) 
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the  Lame '  constants  are  X  and  y,  and  repeated  indices  denote  summation. 


According  to  Snell's  law 
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the  intensities  in  the  x-  and  y-direction  can  be  given  by 
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where  the  asterisk  denotes  the  complex  conjugate,  and  both  ot^  and  are 
either  purely  real  or  purely  imaginary  depending  upon  whether  the  incident 
angle  is  less  than  or  greater  than  the  compressional  or  transverse  Snell's 
law  critical  angle.  To  evaluate  the  wave  amplitude  coefficients  A...F,  we 
have  applied  the  well-known  boundary  conditions  of  continuity  of  displace¬ 
ment  and  pressure  across  the  interfaces  at  x=o  and  x=h,  and  the  shear 
stress  must  vanish  because  a  fluid  medium  on  each  side  of  the  layer  has 
been  assumed.  Mathematically  the  boundary  conditions  can  be  represented 
as 
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and  is  the  Kronecker  delta. 

Solution  of  these  boundary  conditions  results  in  six  simultaneous 
equations,  in  which  the  unknowns  are  the  desired  amplitudes.  Letting 
A...F  »•  Ai...Ag,  the  solution  can  be  put  in  the  form 
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. *  lift  i  7.  u.r  w>.  Jt  . . . 


and 


These  equations  have  been  solved  using  Program  LAYER  in  Appendix  B  for 
the  case  of  sea  ice  bounded  on  one  side  by  water  and  on  the  opposite  side 
by  air.  Although  the  results  of  this  calucation  predict  anomalies  that 
remain  to  be  interpreted,  the  dispersion  of  the  Lamb  waves  can  be  extracted 
from  the  calculation  and  are  shown  in  Figure  30  of  the  text. 
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PROGRAM  RAYQ05(INPUT,OUTPUT,TAPE5=1NPUT,TAPE5»OUTPUT,TAP£99) 

C  RAY  TRACING  PROGRAM, APPLIED  PHYSICS  LABORATORY 

C  RE AO  IN  TABLE  OF  VELOCITIES  FOR  DECREASING  DEPTHS 

DIMENSION  VELOC (508) ,10(500) , DEPTH  <500) ,SHIP<  3) ,OATE (3) 

OATA  TEST1,P1/1. £-6, 3. 141 59265358/ 

CALL  PLOTS 
READ (5, 20 1)  NOFLT 
201  FORMAT (15) 

702  READ (5,896)  SHlP,OATE 

096  FORMAT (3A1Q,3A10) 

DO  3  M=1 ,500 

READ (5,1)  OEPTHCH) ,VELOC (Mi ,ID(M) 

1  FORMAT  (2F20.1,Ili 

IF  <IDCM)-9)3,4, J 

4  II*M-1 
GO  TO  5 

3  CONTINUE 

5  READ (5,6)  OSOURC ,SVEL , RMAX  , DELI  A A , ANGMA A , AN G MIN 

6  FORMAT  <  3F1 0 •  1 , 3F 1 0 • 5  I 

WRITE  (6,60)  SHIP, DATE 

GO  FORMAT  (1H1/1H  34X.20H  RAY  TRACING  PROGRAM/ 1H  31X,27H  APPLIED  PHY  I 
1ICS  LABORATORY/ 1H  30X,3A10/1H  3UX3A1Q/// ) 

WRITE (6, 162)  (DEPTH (Mi , VELOC (Mi ,M=i,II) 

162  FORMAT  <7H  OEPTHaFlQ.l ,10X,10H  VEL QCIT Y  =  F 1 3. 4 > 

WRIT  E  <  6, 1 09 )  SVEL.RMAX  , DELTAA, ANGMAX .ANGMIN 

109  FORMAT  ( 1H  ,  2F  10 . 1 , 3F 1  0. 3/  ) 

READ  IN  THE  SCALE 

READ  (5,7  01)  YMIN,  YMAX,OY,SY,XMIN,XMAX,L)X,SX,DIM,rQ,  VMIN , VN AX , O V , 3 V 
701  FORMAT ( 4F10 . 2/ 6F10 .2/4F1 0 . 2) 

c 

IF (TO.EQ.Q. )  TO  =  DEPTH ( I I ) 

SUM  =0.0 
UO  130  1=1, II 

IF(OEPTH(Li.LT. OSOURC)  GO  TO  130 
IF(DEPTH(L) .GE.TOiGO  TO  119 

SUM=(VELOC(LMVELUC(LPl))/2.*(UEPTH(LU>-OEPTH(L)i  F  SUM 
GO  TO  130 

110  CONTINUE 

IF(OEPTH(L) .ED. TC)  GO  TO  125 
DEL  =  OEPTH(L)-TD 

VIN  *  (VELOC (L-l)-VELQC(L) 1/ (DEPTH  (L-i)-OEPTh (L) I  *  DEL 
SUM  a  SUM  -  (VELOC  (L) -VIN/2. i  ♦  DEL 

125  CONTINUE 
VAV  a  SUM  /  ( TO-OSOURC) 

IF(0 IN. GT • 0 . 5 )  GO  TO  127 
WRITE (6,126)  OSOURC, TO, VAV 

126  FORMAT  ( 1HQ , 25H  AVERAGE  VELOCITY  BETWEEN  ,F8.2»4H  AN0,F8.2, 

X  oH  FEET  IS , FS.2«1UH  FEET/SEC  //28H  ALL  DISTANCES  ARE  IN  FEET 
X  / lHl ) 

GO  TO  135 

127  CONTINUE 

WRITE (6, 128)  OSOURC, TO , VAV 

128  FORMAT  (1HQ,25H  AVERAGE  VELOCITY  BETWEEN  ,F8.2,4H  AND, F8. 2, 10H  MET 

XERS  IS  ,F8.2, 12H  J1ETERS/SEC  //28H  ALL  DISTANCES  ARE  IN  METERS 
X  / lH 1) 

GO  TO  135 
130  CONTINUE 
135  CONTINUE 
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C  PLOT  AXIS 
C 

PLOT  V*  (YMAX-YMIN)/OY 
If <DIM.LT. 0.5)  GO  TO  400 

CALL  AXIS  (  0  »  0  »  0 . 0  ,24H  ,  25, SV, 

13.(1,  VMIN,OV,  10.) 

CALL  AXIS (0. 0 ,PL0TY,13H0EPTH  ( METER) , 13 , SY ,9 0 . 0, YMAX , Q Y,  10.  ) 

CALL  AXIS(0»0,PLGTY, 24HVELQCITY  PROFILE  (M/SEC) ,-25, SV, 

10.0, VMlN,OV,10.) 

l-  CALL  AXIS(SV+1.,PL0TY,13HRANGE  ( METER ) ,-13 »SX  ,  0 . 0 , XMIN  ,QX , 10 , > 

CALL  AXIS(SVH.,PLQTY,  13HOEPTH  (METER) ,13, SY, 90.0, YMAX, OY, 10. > 

CALL  AXIS(SV+l.,a.Q,13H  ,  12, SX, 0 . 0, XMIN, OX, 10 . 1 

GO  TO  410 
400  CONTINUE 

CALL  AXISCO. 0,0.0  ,25H  ,  25, SV, 

10.0, VMIN,0V,10.) 

CALL  AXIS(Q.Q,PLOTY»13  HOE  PTH  (FEET)  , 13 , 3 Y , 9 0 . 0 , YMAX , U Y , 1 0 . ) 

CALL  AXIS(O.OfPLOTY,25HYELOCITY  PROFILE  (FT/SEC) ,-25, SV, 
10.0,VMIN,OV,10.) 

CALL  AXIS (SV+1. , 0.0, 13H  ,  12 ,SX, 0 .0, XMIN, OX, 10 . ) 

CALL  AXIS(SV*l.,PLOTY,13HOEPTH  (FEET)  ,13, SY, 90.0, YMAX, OYilO.) 

CALL  AXISfSV+l.iPLOTY, 13 H RANGE  (FEET)  ,-l 3 , SX , 0 . 0 , XMIN.G X , IQ . ) 

410  CONTINUE 
C 

■  C  WRITE  HEADING 

:  ‘  C 

PL0TX=(VELGC(1)-VMIN)/DV 
!  PLOTY= (-DEPTH ( 1) -YMIN) /oy 

i  CALL  PLOT (PLOTX.PLGTY, 3) 

.00  TOO  IN1=2, II 

PLOTXa  <  VELOC  1IN1)-VMIM)/0V 
PLOTY=(-OEPTH( INI) -YMIN)/ OY 
700  CALL  PLOT (PLOTX,PLOTY, 21 
NOPLT“NOPLT-l 
CALL  PLOT(SVU.,O.C,-3) 

PLQTY*  (YMAX-YMIN)/0Y 
i  PL0TY=PL0TY-.6 

PLOTX=PLOTY-  .4 

CALL  SYM80L<6. 25, PLOTY ,0.14, SHIP  ,0.0,30) 

I  CALL  SYMBOL (6. 25, PLOTX ,0. 14, SATE  ,0.0,30) 

PLOTS Ys(-USOURC- YMIN) /OY 
00  170  1=1,11 

IF  (OEPTH(I)-OSOURC) 170,171,170 
,  171  ISQURC*I 

|  GO  TO  172 

1  170  CONTINUE 

WRITE  (6,202) 

|  ■  202  FORHAT (5X,52H0EPTH  OF  SOURCE  HAS  NOT  BEEN  INCLUDED  IN  YOUR  TABLE  ) 

)  STOP 

1  ■  172  K*ISOURC 

atheta=angmax 

f  :  24  THETAJeATHETA»PI/180. 

I  ANGLEQsASS  (SIN  (THETA J)} 

IF (ANGMIN-ATHETA-TES  ri)9,9,200 

9  WRITE (6,7)  OEPTH(ISQURC)»ATHETA 

?  FORMAT  (14HCS0URCE  0£PTH*F8.2,  3X.13HS0URCE  ANGLE=F«. 2, 8H  CEG 

1REES//4X5HPQINT5X  7HRANGE  3X10HT IME  (SEC)  6X10HOEPTH  ///) 

‘  CALL  PLOT (0.0,PLQTSY,3) 

10  R*0. 

,  IP*l 

!  161  T*0« 

IPOINT=l 

11  RL*COS  ( THETA J)/ VELOC ( IS OURC) 

/  12  IF  (IP-1) 13,14,13 
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14  H0LD*A6S  (THETAJ) 

IF  (H0LD-(i.QE'6)l 330,821,821 
821  IF  (THETAJ) 18,330,500 
500  I*ISOURC  -1 
GO  TO  13 

330  J*ISQURC-1 
N=ISOURC 

AHIGH=(VElOC(J)-VELUC(N)  )  /  A03  l DEPTH ( J ) - OE PTH l N) ) 

J»ISOURC*l 

ALOW  =<VELOC(  J)-VELQC(N)>/A8S  ( DEPTH (J ) -DEPTH <N)  ) 

PLOTX=RMAX/DX 

332  IF  (AHIGH)333, 331,150 

333  l=ISOURC-l 
IZ£RO=l 

GO  TO  13 

150  IF  (AL0W)151, 152,153 

151  IZER0=3 
I*ISOURC*l 
GO  TO  13 

153  IZEROsS 
GO  TO  22 

152  WRITE (6,  334) 

334  FORMAT  (55H0  RAT  TRAVELS  HORIZONTALLY  TO  NEXT  PROFILE 

I///) 

IZER0=3 

CALL  PLOT (PLOTX.PLOTSY ,2) 

GO  TO  22 

331  WRITE (6,334) 

CALL  PLOT (PLOTX.PLOTSY, 2) 

IF  (ALOW  1  335,153, 153 

335  WRIT E  t  6, 336) 

336  FORMAT  (64HOZERO  GRADIENT  ABOVE , NEGATIVE  GRADIENT  BELOW. RAY  ALSO  3 
1ENOS  OQWN///> 

GALL  PLOT(Q.O,PLOTSY,3) 

GO  TO  151 

16  I»ISOURC*l 

13  IF  (1)17,18.17 

17  IF  (11-1)19,20,20 

18  WRIT  £  ( 6 , 21> 

21  FORMAT  (29H  RAY  GOES  ABOVE  MINIMUM  DEPTH//) 

GO  TO  22 

19  WRITE (6,23) 

23  FORMAT  (29H  RAY  GOES  BELOW  MAXIMUM  DEPTH//) 

22  IP=1 

ATHET  A*A YHET  A-OELTAA 
ANGLEBcAQS  (SIN  (THETAJ)) 

K»ISOURC 
IZER0*3 
GO  TO  24 

20  VL*VELOC ( I) *RL 

£5  IF(VL.LE.l.)  ANGLEA*SQRT(t.-VL**2» 

26  Asi'VELOC  (I)-VELOC(K)  )/ABS  (DEPTH  (I  J-OEPTH  (K) ) 

27  IF  ( A ) 28 , 29, 28 

28  IF  (VL-1. >30,30, 31 

31  ANGLE A=G. 

IN*i 

GO  TO  32 
30  IN*3 

32  OELTAR=(ANGLE9-ANGLEA) /A/RL 

33  UELTAT«AL0G ( ( 1 .♦ ANGLES )/ ( 1.- ANGLES )* ( 1 .-ANGLE  A)/ (l.*ANGLlA))/2,/A 
IF  (IN-1) 37,34,37 

34  D=tl,/RL-VELOC(K))/(VtLOC(I)-VELOC(K)  ) 

0*0*A9S  (CEPT  H ( I ) -DEPT  H (K) ) 

ISAVE*K 
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36  IF  (1-10  33,200,39 
39  l)U=0£PTH ( <) +0 

INEXT=K-1 
GO  TO  41 

36  DO*Ci£Pl  H  ( K)  -0 
INEXT=K+1 

41  R=R* A8S  (QELTAR) 

T*T*ABS  (OELTAT) 

PLOTX=R/OX 
PLOTY=(-00-YMIN)/UY 
CALL  PLOT  (PlCTX.PlOTY, 2) 

IP  (ATHET A. £Q.- 6.. ANU.PLOTY.LT. -2. )  CALL  SYMBO L CPLOTX , PLO T Y , . 0 15 , 
114,0., -2) 

WRIT  £(6,42)  IPOINT,R,T,OQ 

42  FORMAT (5XI3,5XF7.1,5XF8.5»5XF8.2) 

43  IF  ( X  N-l >44,45  ,44 

45  R=R*ABS  (UELTAR) 

T*T+A9S  (OELTAT) 

IPOI NT=I POINT  *1 
I*ISAVE 

IN»2 

GO  TO  46 

37  R»R*ABS  (DELTAR) 

T*T*A6S  (OELTAT) 

46  WRITE (6,42)  IPOINT , R , T , DEPTH C I ) 

PLOTX-R/OX 

PLOTY=(-DEPTH(I)-YMIN)/OY 
CALL  PL0T(PLCJTX,PlQTY,2) 

IF <ATHETA.EU. -6. .ANU.PLOTY.LT. -2.)  CALL  SYMBOL CPLOTX , PLO I Y ,. 0 15 , 
114, Oo, -2) 

GO  TO  43 

29  OELTARxAOS  (OEPTH(I)-OEPTHiK) ) *VL/ANGLEA 
46  OELTAT  =ABS  IOEPTH (I) -QEPT H C K ) )/tf 6L0C <1 >/ ANGLE A 
IN*3 

GO  TO  37 

44  IF  (IN-2)  49,50  ,49 

49  IF  CI-Kl 53,200,52 

53  K=I 
1*1-1 

GO  TO  54 
52  K*I 
1*1*1 

54  IPOI NT*I POINT + 1 
ANGLEdsANGLEA 
GO  TO  55 

50  K* ISA  YE 
IPOINT*IPOINT*1 
ANGLEi3*ASS  (ANGLES) 

I-INEXT 

55  IF  <RMAX-R)56,57,57 

56  IF  (I2ERO-l)137,136,137 
137  WRITE  <6, 611) 

611  FORMAT  1 26H  RANGE  EXCEEOS  MAXIMUM  RANGE///) 

GO  TO  22 

57  IP*2 

GO  TO  12 

136  IF  (ALOWI 139, 152,153 

139  WRITE  (6 , 140 ) 

140  FORMAT  ( 54M  NEGATIVE  GRADIENT  BELOW  AND  ABOVE. RAY  ALSO  SENDS  DOWN/ 
l//> 

I*ISOURC  *1 

IZERO=3 

IP*2 
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ANGEE=THETAJ* 57. 295/8 
ANGLEB=A3S  (SIN  (THtTA J) ) 

KaJSQURC 

WRITE (6, 7)  OEPTHCISOORC) .ANGLE 

R*Q. 

GO  TO  161 

200  FtNAL  =  SXHO.  0 

CALL  PLOT (FINAL, Q. 0,-3) 

NOPLT  =  NOPLT-1 
IF (NOPLT . NE » 0 )  GO  TO  71)2 
CALL  PLOT (0 • . 0 . . 339) 

STOP 

ENO 


PROGRAM  SSCTR (INPUT, OUTPUT, TAPE5= INPUT, TAPE6»OUTPUT) 

C  THIS  PROGRAM  CALCULATES  THE  SCATTERING  CROSS  SECTION  TIMES  THE 

c  frequency  squared.  the  value  is  calculated  for  a  shear  plane 

C  WAVE  IN  A  SOLID  STRIKING  A  FLUID  FILLED  SPHERE.  THE  CASE  OF  A 

C  COMPRESS IQNAL  WAVE  STRIKING  THE  SPHERE  IS  ALSO  CALCULATED.  THE 
C  VALUE  IS  CALCULATED  AS  4  FUNCTION  UF  THE  U IMENS IONLESS  PARAMETER 

C  KA  WHERE  K=2*PI*FR£Q/ VELOCITY  AND  A  IS  THE  SPHERE  RAOIUS. 

C  THE  CALCULATION  IS  MAO E  IN  THE  C . C.S. SYSTEM.  THE  SCATTERING 
C  CROSS  SECTION  HAS  UNITS  OF  0B-CH**2.  CS1  IS  THE  SHEAR  VELOCITY 

C  AND  CL1  IS  THE  COMPRESSIDNAL  VELOCITY  IH  THE  SOLIO.  CL2  IS  THE 

C  COMPRESS IONAL  VELOCITY  IN  THE  SCATTERER.  RHO  IS  THE  RATIO  OF  THE 

C  DENSITY  OF  THE  SPHERE  TO  THE  DENSITY  OF  THE  SOLID. 

READ (5,49)  CS 1 , CL1 , CL2  ,RHO 
49  FORMAT (3FS.Q,F6.3) 

WRITE (6, 51) 

SI  FORMAT (42H  KA  GAMMA, SHEAR  GAMMA, COMP) 

1  READ (5 ,5  0 )  AKA 
50  FORMAT (F?. 2) 

XF(EOF,5)  7,8 
8  CONTINUE 

PI*3. 141592654 

C  FIRST  THE  SHEAR  SCATTERING  IS  CALCULATED. 

BKS1»AKA 

3KL1=3KS1*CS1/CL1 
3KL2*BKS1*CS1/CL2 
SIH*CS1/ (2. *P I ) 

C  IF  CIN=0.  Then  A  SHEAR  HAVE  IS  INCIDENT. 

C I H3  O  . 

CALL  COEF ( 8KS 1, 3KL 1, 3KL2, SIM, Cl M, RHO, GAMAS) 

GAKAS»GAMAS*3KS1 

C  GAMAS  IS  THE  SHEAR  SCATTERING  CROSS  SECT ION*F REQ**2 . 

C  NOW  THE  COMPRESS IONAL  SCATTERING  IS  CALCULATED. 

3KL1“AKA 

aKSlaQKLl*CLl/CSl 
3KL2*BKL1¥CL1/ CL 2 

C  IF  SIM=Q.  THEN  A  COMPRESSIDNAL  WAVE  IS  INCIDENT. 

SIHaQ. 

CIM*CL1/ (2.*P I ) 

CALL  COEF ( QK51,QKLi,0KL2,SIM,CIM, RHO, GAM4C) 

GAMAC*GAMAC*8KL1 

C  GAMAC  IS  THE  COMPRESSIONAL  SCATTERING  CROSS  SECTION*FREQ**2. 

WRITE (6, 55)  AKA, GAMAS, GAMAC 
55  rORNAT(2XF7.2,Z(5XE13.4>) 

GO  TO  1 
7  CONTINUE 
ENO 
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SUBROUTINE  COEF ( AKS1 , AKL1, AKL2, SI , Cl , RHA , GAMM A) 

OIMENSION  AJ1 (400) ,  AJ2(400>,  AJ3(4Q0),  Y11430),  Y2<400) 

THIS  SOLOES  THE  SOLING  ARY  CONDITIONS  FOR  THE  AMPLITUOES  OF  THE 
SCATTERED  HAVES.  THE  AMPLITUDES  ARE  USED  TO  CALCULATE  THE 
SCATTERING  CROSS  SECTION. 

CALL  SPBES(AKS1,N1,AJ1) 

CALL  SPBES(AKL1,N2,AJ2) 

CALL  SPBES(AKL2,N3,AJ3) 

PI*3. 141592654 
IF(N1-N2)  10,11,11 

10  N=N2 

GO  TO  12 

11  N»N1 

12  IF (N-N3)  13,14,14 

13  NM*N3 

GO  TO  15 

14  NH*N 

NH  IS  THE  UPPER  LIMIT  OF  THE  SUM  IN  THE  SCATTERING  CROSS  SECTION. 

15  CALL  SNHN(AKS1,NH,Y1) 

CALL  SNMN(AKL1,NM,Y2) 

G AMHAz  0 . 

00  100  M-l , NM 
F*FLOAT ( M) -1 . 

ZERQ*0. 

011sF*AJ2(M)-AKl1*AJ2 (M+1) 
l)12=ZEK0-F*(FU.  )*AJ1  (M) 

RAT»F* (l.-F) +AKSl**2/2. 

021=RAT*AJ2(M) -2  .  * AKL 1 *A J2 (M+ 1 ) 

022=*F*(F*1.)* < <F-1.)*AJ1(M)-AKS1*AJ1(M+1) ) 
031*(F-1.)*AJ2<M)-AKL1*AJ2(M+1) 

ROT»l.-F**2+AKSl**2/2. 

Q32*R0T*AJ1(M)-AKS1*AJ1(N+1) 

CI11SF*Y2(M)-AKL1*Y2(M+1) 

CI12=F*(F+1.)*Y1(M) 

C13*F+AJ3 (M) -AKL2* AU3 ( M+l) 

CI21zRAr*Y2(M>-2.*AKLl*Y2(MU) 

CI22*F*(F+1.)  *(<F-1.  I  ♦  VI  (M)-AKS1*Y1<M+1) ) 
C23=RHA*AKSl**2*AJ3(M>/2. 

CI31UF-1.)  *Y2  (M)  -AKL1*Y2 (M+1) 

0132*  R0T*Y1(M)-AKS1*Y1(M+1) 

Z1»(011*CI<012*SI) 

Z2»(021*CI+022*5I) 

Z3*(03i*CI+032*SI) 

ARM*022*C13*Z3+C23*032*Z1-C23*Z3*012-D32*Z2*C13 
AIM*C13*CI22*Z3+C23+CI32*Z1+C23*Z3*CI12-CI32*Z2*C13 
aRM*031*Z2*C13+C2  3*Z3*QU-C13*Z3*D21-C2  3*Zl*D31 
OIH*C13*Z2*CI31+C23*Z3*CIil-Ci3*Z3*CI21-C23*Zl*GI31 
0RH«C23*(012*031*CI12*CI 31-D11*Q32+CI11*CI32)  +C13*  (021*0  J2-CI21* 
10132-022*031+ CI22*CI31) 

OIM»C2  3*  (Cl 3 1*012-0112*0 3 1-CIU  +  Q32-C  132*0 11) +Cl3*<D32*CIZi  +  D21* 
1C 1 32 -C 122 *03 1-02 2* Cl  31) 

P*ARM/ORM 
Q*BRM/ORM 
R»OI M/URM 
SsAIM/ORH 
T*8IM/0RM 
UENz 1 . +R**2 

C  AM  ANO  BN  ARE  THE  SQUARES  OF  THE  MAGNITUDE  OF  THE  SCATTERED  HAVES 
AM*<P**2+S**2)/0EN 
BH»<  Q**2  +  T**2 ) /QEN 

GP»( 2.*F  +  1. )  *  4.*PI*<  AM/ AKL1+F*  (F+1.)*BM/AKS1> 
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.‘ft’.-;  .«  ..  .  I'.'i.i  ... 


oo  nnn  o.o  o 


GAMNA=GAMMA  *GP 
IOC  CONTINUE 

G A MM A *0. 686* GAMMA 

RETURN 

ENG 


SUBROUTINE  SNMN (ZT «NN«  Y) 

DIMENSION  Y1400) 

THIS  CALCULATES  THE  SPHERICAL  NEUMANN  FUNCTION (Y < I ) >  OF  4RG  ZT  AND 
UP  TO  ORDER  NN.  THE  VALUES  ARE  CALCULATED  BY  UPWARQ  RECURSION 
RELATIONS. 

ZER0*Q. 

DO  1  I«lc400 
1  Y (I) *0. 

Y ( 1)  = ZERO-COS  ( ZT ) /ZT 

Y (2)  *ZERO-COS (ZT)/ZT**2-SIN(ZT)/ZT 

NT»NN*2 

DO  2  1=3, NT 

AJ=FLOAT ( I) -2. 

Til) =( <2.*AJ+1,)/ZT>*Y <1-1 >-Y (1-2) 

2  CONTINUE 
RETURN 
END 


SUBROUTINE  SP8ES (ZT, NT , A J; 

DIMENSION  R ( 4  0  C  >  >  RJ<  400),  A  J (400) 

THIS  CALCULATES  THE  SPHERICAL  BESSELS  FUNCTIONSCAJ(I) )  OF  ARG 
ZT  UP  TO  ORUER  NT, 

THIS  ROUTINE  IS  VALID  FUR  ARGUMENTS  AS  LOU  AS  0.05. 

BT=ZT 
NOIH*400 
NZT» IFIX (ZT ) 

DO  100  1  =  1, NO  I M 
100  AJ(U=0. 

THIS  SECTION  SETS  THE  UPPER  LIMIT  OF  NT  FOR  THE  GIVEN  ARG  SUCH 
THAT  ALL  AJ*S  OF  HIGHER  ORDER  CAN  BE  SET  TO  ZERO. 

N=NZT*10 

OO  1  I=N,1B5Q,5 

J  =  I 

FI*FLOAT (I) 

StCA=ZT/ (FI *. 5 ) 

TACA=SQRT  U.-SECA**2J 
GOSA»l./SECA 
SINA=SQRT (l.*COSA**2) 

ALP=ALGG (COSA tSQRT (COSA**2-l.  ) ) 

0£LT=EXP( (FI+.5)*(TAGA-AlP)>/ ( 2. * < FI+ .5 ) *SQRT (SINA) ) 

IF  (OELT-l.E-27)  5,5,1 
1  CONTINUE 
GO  TO  12 
5  N»J 

FN=FL0AT(N) 

NP»N*1 

NT*N+2 
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NU*N*lO 
EXPR=EXP(1.1 
00  7  I=NU ,1890,5 
K*I 

FI*FLQAT ( I) 

0*0.434294 

Al=2.* (FI-FN+1 .J*0*AL0G(ZT) 
A2=<2.*FI-2.*FN+1.)*0*ALOG<EXPR) 

A3*lFN+2. 5)  *0*  ALOG  <FN*  2. )  ♦  £FN-1 . 5 ) *0* ALOG (FN- 2. ) 
A4*< FI ♦3.51*0* ALOG(FI+  3. ) 
A5=<FI-.5S*0*ALOG(FI-1.J 
E=A1+A2+A3-A4-A5 
IF  (£>1Q.  )  9,9,7 
7  CONTINUE 

THIS  SECTION  CALCULATES  THE  BESSEL  FUNCTIONS. 

GO  TO  12 
9  NU*K 
j=nu 

R<NU+ll=0. 

00  2  1=1, NU 
K*I 

SJ=FLOAT<J! 

R(  J)  -ZT/ <l.+2.*SJ“ZT*RU*in 
IF(R(J>-1.)  2,2,3 

2  J*J-1 

GO  TO  200 

3  IF<J-2>  200,2000 14 
200  RJ<NU*1)=R(NU» 

RJ(NU) =1. 

J— NU 

NUP=NU-1 

00  210  1=1, NUP 

J*J“1 

SJ*FLOAT  CJ) 

RJ(J)=ll.t2.*SJ)/ZT*RJ (Jfi>-RJ(J*2> 

210  CONTINUE 

ALPH=(RJ<i)-ZT*RJt2)  1  *  COS  ( 3mZT*  S  IN  ( BT  >  *  R  J  <  1  > 
00  220  1=1, NT 
220  AJ(IJ=RJ(I)/ALPH 
333  CONTINUE 
RETURN 

14  RJIJ*1>=R(J) 

LAH2*J*2 

IF<LAN2-iNF)  15,44,44 

15  RU (J ) =1. 

«i*J-l 

K*K»1 

00  3  I*K , NU 
SJ*  FLOAT (J) 

RJCJJ*ti.*2.*SJ>/ZT*RJ  £ J+l J-R  J  C J  +  2 ) 

00**4*  I  =  L  AM2 ,  NT 
4  RJ (I J  *RJ ( 1-1 S  *R ( I- 1) 

ALPH* (RJ(1I-ZT*RJ(2£) *COS { 8T) *ZT* SIN t BT) *RJ < 1 ) 
NUH*N0IM-1 
IF  <NT~NO  HI  16,16,17 
17  NT*NOH 

16  00  6  1*1, NT 

6  AJ(I)=RJ(I)/ALPH 
105  CONTINUE 
RETURN 

44  HRITE (6, 45) 

45  FORNAT  (25H  LA M BOA  +  2.  GE  .(£♦  1  END  XEO) 

GO  TO  1Q5 


non  o  on  c  o  oooo  ooooooo  o  rs  n  o  o 


12  WRITE (6,13)  ZT 

13  FORMAT (26H  ARG  BES  FTN  TOO  LARGE  ZT=,E20 .8 ,5X14H  ENO  EXECUTION) 
STOP 

ENO 


PROGRAM  LAYER  < INPUT ,  OUTPUT  » TAPES  =  INPUT  *TAPE6=QUTPUT ) 

THIS  PROGRAM  SOLVES  THE  PROBLEM  OF  AN  ELASTIC  LAYER. 

BOUNDED  ON  EACH  SIOE  dY  OIFFERENT  IDEAL  FLUIDS. 

IMPLICIT  COMPLEX (A, U,F,D,P»S! ,REAL(K> 

IMPLICIT  COMPLEX  NOTATION  MAY  BE  VALID  FOR  THE  U/W  COC  640Q  ONLY. 
OIFFERENT  COMPUTERS  MAY  REQUIRE  ANOTHER  NOTATION. 

COMPLEX  CC,0O,EE*GAB,CSINE»Q,J 
01  MENS  ION  F (6 )  »  KK  ( b) 

NECESSARY  FUNCTION  STATEMENTS  FOLLOW. 

CAB(A,B)=A*C0NJG(8) 

ARG< A,B) =0. 5* ( AfCDNJG<di  ) 

CSINE<A,a)=CSINtARG(A,ai )/ARG(A,9) 

J*CMPLX ( 0 . i 1. ) 

DATA  INPUT  REQUIREMENTS  TO  THE  PROGRAM  FOLLOW- 

NNN  IS  THE  NUMBER  OF  INDIVIDUAL  SETS  TO  BE  CO NSIOE REO , CO ,CL , CT , AN ) 
CF  ARE  THE  ACOUSTIC  WAVS  VELOCITIES  OF  THE  INCIDENT  WAVE,  THE 
LONGITUDINAL  AND  TRANSVERSE  WA YE  IN  THE  PLATE  AND  THE  VELOCITY  OF 
THE  TRANSMITTED  WAVE  RESPECTI VEL Y .  RHO.RHOO  AND  RHOF  ARE  THE 
DENSITIES  OF  THE  PLATE,  INCIDENT  MEDIUM  AND  FINAL  MEDIUM 
RESPECTIVELY  IN  GRAMS  PER  CUBIC  CENTIMETER. 

RE AO ( 5 ,99)  NNN 
99  FORMA T (1 2 ) 

NNX*1 

10  0  READ  15,101) CO , CL ,CT , CF ,RHO, RHOO , RHOF 

101  FORMAT  <4F5»  0 »  3F7.4I 

KH  IS  H*P I*  (T  HICKNESS/ WAVELENGTH  IN  THE  INCIDENT  MEDIUM) ,CHI0  IS 
THE  INITIAL  INPUT  ANGLE, OCHI  IS  THE  INCREMENTAL  ANGLE  OF  CHANGE, 
CHIP  IS  THE  FINAL  ANGLE  TO  BE  CONSIDERED  AND  X  IS  THE  DEPTH/WAVE¬ 
LENGTH  IN  THE  INCIDENT  MEDIUM. 

102  REAO (5,401) KH ,CHI0,DCHI, CHIP, X 
401  FORMAT (5F8.4) 

PI»CMPLX (3.141592654,0.) 

R»RHO/RHOQ 

R1=RHOO/RHOF 

DETERMINATION  OF  THE  ACOUSTIC  RUNNING  VECTORS  FOLLOW- 

KLH*(CO/CLl*KH 

KTH*(C0/CT)*KH 

KFH=*  (CO/CF)  *KH 

LL  DETERMINES  THE  NUMBER  OF  ANGULAR  INCREMENTS  TO  BE  CONSIDERED 
LL*IFIX((CHIF-CHIQ)/OCHI) 

ALTHOUGH  PI  IS  SPECIFIED  COMPLEX,  ONLY  THE  REAL  COMPONENT 
CONTRIBUTES— STATEMENTS  CHANGE  DEGREES  TO  RADIANS 
CHIO*CHIC*PIr 180. 

OCHI*OCHI*PI/18Q. 

FOLLOWING  STATEMENT  SETS  INITIAL  ANGULAR  PARAMETER. 

CHI*GHI0 
00  121  IL»1,LL 
IF(CHI)200,201,200 
200  IF(IL-l) 202,202,201 
2C1  CHI«CHI+OCHI 

FROM  SNELLFS  LAW 

COMPLEXITY  OF  THE  FOLLOWING  STATEMENT  ORIGINATES  FROM  LACK  OF 
INTESTINAL  FORTITUDE  TO  USE  BET A*KH*CS IN < CHI I . 

202  BETA»KH*CSIN<CMPLX<CHI ,0. ) ) 


C  DIRECTION  COSINE  DETERMINATION  FOLLOWS- 
AOH=CMPLX(KH*COS(CHI) , 0. ) 

C  CONJUGATE  CONSIDERATIONS  NECESSARY  TO  ASSURE  A  OECAYING  WAVE  IN 

C  THE  MEDIUM  FOR  ANGLES  OF  INCIDENCE  BEYOND  THE  CRITICAL  ANGLE. 

ALH=C0NJG<CSQRT(KlH«*2“3£TA**2) > 

ATHaC0NJ&(CSQRT<KTH**2-3£TA**2) ) 

AFH*CQNJG (CSQRT ( KFH**2-3E  T  A**2  > ) 

IFCCF-C0)lQ8,107,loa 

107  AH=(1.,0.) 

GO  TO  109 

108  AH=A£1H/AFH 

C  THE  SPECIFICATION  OF  AFH  AND  R1  ASSUMES  DIFFERENT  MEDIA  ON  EACH 

C  SlOt  OF  THt  PLATE  IF  DESIRED. 

109  S12=-2.*ALH*8ETA/(ATri**2-BETA**2) 

S2Z*ALH/A0H 

S24aSETA/A0H 

S32»R*  (2.*ttET  A**2/KTH**2-1.) 

S34a-R*2.*BET  A'*ATH/KTH**2 
S1*(S12*S24-S22) *AH 
Glal.-CCOS(ATH)*CCOS (ALri) 

G2=CSIN(ATHI*CSIN(ALH) 

G3*CC0S( ATM J*CSIN( ALH) 

G4aCCQS(ALH)*CSIN(ATH) 

A1=4.*(2.*K1*S12*SJ2*S34*G1-(R1*  t  <  S12*S  34)  «  *  2  «-S32**2  > 
1-S1**2)*G2+J*  <S12*S24-S22)  *<4H-R1)  ♦ ( -S12*S 34* G3+S 32*G 4 > ) 

A2=2.*  <R1*S12*S34*G1-(R1*S32*-Si)  *G2+ J* ( R 1*S32 +S1 ) *G4 
1-J*R1*S12*S34*G3) 

A3»2.*  C-R1*S12*S34*G1  +  <R1*S32-S1 1  *  G2  *•  J*  l  R1*S32-S1 )  *G4 
1-J*R1*S12*S34*G3) 

A4=2.*S12*(R1*S32*G1-(R1*S12*S34-S1)  *G2+J*(RI*S12*S34-S1)*G3 
l-J*Ri*S32*G41 

A5=2.*S12*(R1*S32*G1-(R1*S12*S34*S1>*G2-J*(R1*S12*S34*S1)*G3 

1*J*R1*S32*G4) 

A6--8.*J*R1*AH*S1* <S12*S34*CSIN(ALH)-S32*CSIN (ATH) ) 
Da4.*<2.*Rl*S12*S32*S34*Gl-<Rl*<  < S  12*S34 > **2» S32**2)  +  S 1* *2) *G2 
1*J*<S12*S24-S22)*(AH*R1)*(-S12*S34*G3*S32*G4) ) 

F  Cl) *A I/O 
FC2>a-2.*A2/U 
FC3)*2.*A3/0 
FC4) »-2.*A4/0 
F  C  5) *2 .* A  5/D 
FC6)a-A6*CEXP( J*AFH)/0 
00  203  IJ=1 i 6 
203  KK  (I J ) =CABS  <  F (IJ) ) **2 

C  THE  FOLLOWING  DETERMINES  THE  SNELLAS  LAW  CRITICAL  ANGLE  FOR  THE 

C  COMPRESS  ION A L  ANO  TRANSVERSE  ACOUSTIC  WAVES  RELATIVE  TO  ANGLES  OF 

C  INCIDENCE  IN  DEGREES. 

critl*asin<co/cl)* lao./Pi 

CRITT-ASIN!C0/CT)*180,/PI 

THETAsCHI-UdQ./PI 

C  THE  FOLLOWING  CHECKS  TO  DETERMINE  IF  ENERGY  IS  CQNSERVEO. 

OH EK* REAL  CKK  C 1 >  +  CAFH/ A  OH) *KK  C  b) /R1 ) 

C  IT  IS  NOT  OBVIOUS  THAT  UNITS  ARE  CONSERVED  IN  THE  FOLLOWING 

C  STATEMENTS--HOHEVER,  RECALL  THAT  X  IS  THE  (DEPTH  INTO  THE  PLATE) 

C  OIVIOEO  BY  (WAVELENGTH  IN  THE  INCIUENT  MEOIUM). 
ARGla-J*2.*PI*ALH*X/KH 
ARG2aARGl*ATH/ALH 

C  THE  FOLLOWING  Ub TtRMINES  THE  COMPRESS IONAL  AnD  TRANSVERSE 
C  DISPLACEMENT  POTENTIAL  REQUIREMENTS  OF  THE  PROBLEM. 
PHIN3F(2)*CEXPCARG1) 

PHIP*F(3)*CEXP<-ARG1> 

PSIN»F(4)«CEXP(ARG2) 

PSIP*F (5) *CtXP (-ARG2) 
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PHI = PH IN* PH IP 
PSI*PSIN*PSIP 
PHN*PHIP-PHIN 
PSN*PSIP-PS1N 

THE  FOLLOWING  DETERMINES  THE  ENERGY  PER  UNIT  AREA  PER  UNIT  TIME 
( INTENSITY)  FLOWING  IN  A  DIRECTION  PERPENDICULAR  TO  THE  PLATE 
SURFACE  (MULTIPLIED  BY  THE  SOUND  VELOCITY  IN  THE  INCIQENI  MEDIUM) 
DIVIDED  BY  THE  (DENSITY  OF  THE  PLATE  MULTIPLIED  BY  THE  FOURTH 
POWER  OF  THE  CIRCULAR  FREQUENCY). 

GX*( l./KH)*  (  ARG (ALH,  AL  H)  *  ( C  AB  (  PH  IN,  PH  IN) -C  AB  (  PHIP ,  PHI  P) ) 
1*ARG(ALH,-ALH ) * ( CA 3 < PH  IN, PHI P ) -C AB (PHIP,PHIN>  ) 

2*ARG(ATH,ATHt*  (CAD  (PS  IN,  PS  IN)  -CA  8  (  PS  IP  ,PS  IP  >  ) 

3*  ARG (ATH,-ATH) *(CAd(PSIN,PSIP>-CA8 (PSIP.PSIN) ) > 

Bl=l.+2.*(3tT  A  **  2  +CAQ  (ALH,  ALH)-KLH**2  )/(KTH**2  ) 
92*2.*(CAa(ATH,ATH)-Q.5*(ATH**2  -BETA**?  >)/(KTH**2  ) 
Q3a(2.*B£TA**2  -KLH**2  *0.5*<TH**2  )/(KTH**2  ) 

84=(2.*BETA**2  -O.S*KTH**2  )/(KTH**2  ) 

THE  FOLLOWING  UETERHINES  THE  ENERGY  PER  UNIT  AREA  PER  UNIT  TIME 
(INTENSITY)  FLOWING  IN  A  DIRECTION  PARALLEL  TO  THE  PLATE  SURFACE 
MULTIPLIED  BY  THE  SOUND  VELOCITY  IN  THE  INCIDENT  MEOIUM  OIVIDED 
BY  THE  FOURTH  POWER  OF  THE  CIRCULAR  FREQUENCY. 

GY»< l./KH)* (OETA*Bl*CAB(PHI,PHI) *BET A*B2*CAB ( PSI > PSI ) 
1-(4.*BETA/KTH**2.) * ( C A B l ALH , A LH) * (CA3 (PHIP, PM  IN) *CAB ( PHI N , PHIP ) ) 
2*CAd(ATH,ATH)  «  (CAB  (PS I P ,  PSI  N>  +CA3  ( PSI N , PSIP )  )  ) 

3*B3* ( C  AB ( PS  N i PHI ) * ATH+ C AB ( PHI , PSN) *C0NJG(ATH) ) 

4-B4*(CA6(PHN,  P5I)*ALH+CAB.(PSI ,  PHN)  *C0NJG  (ALHI  )  ) 

IF  (CHI -REAL  US  IN  (CO/CL)  )  )  3  30,  30  0,301 

300  BN*1 

BPsC'SINE  (ALH,  ALH) 

PN*1 

PP=*CSINE(ATH,  ATH) 

GO  TO  304 

301  BN=CSINE(ALH,-ALH) 

BHsl 

IF(CHI-KEALIASIN(CQ/CT)>) 302,302,303 

302  PN-1 

PPaCSINE(ATH, ATH) 

GO  TO  304 

30 3  PN*CSIN£(ATH, -ATH) 

PP  =  1 

304  CONTINUE 

THE  FOLLOWING  SPECIFIES  AMPLITUDES  AND  ALGEBRAIC  TERMS  NECESSARY 
FOR  THE  DETERMINATION  OF  THE  TOTAL  ENERGY  BEING  PROPAGATED 
PARALLEL  TO  THE  PLATE  BOUNDARIES  AS  A  FUNCTION  OF  THE  INCIDENT 
ANGLE, I. SIGMA  IS  THE  INTEGRATED  INTENSITY  ACROSS  THE  PLATE. 
AA*F(1) 

aa*F(2) 

CC*F(3) 

00*F(4) 

EE*F ( 5 ) 

FF«F(6) 


Pl*(CA3(QB»dB)*CEXP(-J*ARG(ALH,-ALH) ) *C A  8 <CC, CC ) *CEXP ( J* ARG ( ALH, 
i— ALH)  )  }*QN*iCAB(BU»CC)  *CEXP  (-  J*  A  RG  (ALH ,  ALH  i  i  CA3  ( CC ,  33 )  *  Cc  XP  (  J* 
2ARG( ALH, ALH) )  )*BP 

P2*(CA0(OO,DO) *CEXP(-J*ARG ( ATH,- ATH) ) *CAB (EE , EE ) *CEXP ( J* 1RG< ATM, 
1-ATH) ) )*PNMCA6(0D,EE) *CEXP (- J* ARG (ATH, A TH ) ) *CAB(EE,00)*CEXP( J* 
2ARG( ATH, ATH) )  ) *PP 

P3*(CAB(86,CC) *CEXP(-J*ARG(ALH,ALH) ) *CAB (CC,OB)*CEXP(J*ARG(ALH, 
1ALHI ) > *BP 

P43(CAB(OD,EE (*GEXP(-U*ARG(ATH,ATH) ) *CAB (EE,O0)*CEXPCU*ARG (ATH, 
1ATH) ) ) *PP 

P5*(CA8 ( EE, 88 ) *CEXP( J* ARG l ATM , AL  H ) )-CAB(QQ,CC) *CEXP(- J*ARG (ATH 
1 , ALH) ) )*ATH*CSINE(ATH,  ALH) 

P6=(CAa(Q0,8BI*CEXP(J* ARG(-ATH,ALH))-CA8(EE,CC)*CEXP(-J**RG(-ATH 
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1 , ALH )  )>*ATH*CS1N£(-ATH,ALH) 

P7=<CAB(CCf UO)*CEXP( J*ARG< ALH, ATM ) ) -CAB < BB ,EE > ♦CEXP (- J*ARG (ALH 

I,  ATM) ) )*ATH*CSIN£(~ATH,ALH) 

Pa-<CAB(CC,££)-»CEXP« J*ARG(ALH,-ATH) > -CAB <B6 , QD) *CEXP ( - J* ARG ( ALH 
l,-7.rH)  )  ) *ALH*CSINE  tALH,ATH» 

SXGMA=B£T A* (dl*Pl+B2*P2-A.*(CAB (ALH, ALH)*P3+CAB( ATH, ATH) *P4> ) 
1*2.*B3*  (ARC  (P5,P5)-ARG  ( P6  ,  P6  »l -2  .  *  B4*  ( ARG  (  P6,  P6>  ARG  (  P  7,  a  7)  > 

SPEC  IF  rNG  REAL  IN  THE  FOLLOWING  ONLY  SERVES  TO  ACCQMOQ AT i  MINIMUM 
OUTPUT  telOS  AS  THE  IMAGINARY  COMPONENT  OF  GX.GY,  AND  SIGMA  MUST 
8E  ZF.Ri  TNCE  THESE  PARAMETERS  ARE  REAL. 

305  KX=REAL  X) 

KY=R£AL(GY> 

TENGY=REAL(SIGMA) 

VR=C0*KH/REAl (BETA) 

IF (IL-1) 116,116, 11 

116  WRIT £(6,204) 

204  FORMAT ( 1H1 ) 

WRITE  (6, 205)  CO  ,  C'. ,  CT  ,  RHO ,  KH,CRITL ,  CRITT 

205  FORMAT  <  3X , 5HC  0  =  , F5. 0 , 3X , 5HC  L  =  , F5. 0 , 3X , SHOT  =  ,F5.0, 

J. 4X  ,6HRHO  -  ,F4.2,3X,6HKQH  =  ,  F4. 1 , 3X  ,  13HPHIL  (  CRI T  )'  -  ,F6.3, 
23X,13HPHIT (CR1T)  =  ,F6.3,/> 

WRITE  (6,117) 

117  F0RMAT(5X,5HTHETA,5X,2HAA,6X,2HBB,6X,2HCC,6X,2H0D,6X,2HEE, 
lt)X,2HFr ,3X,5HCH£CK,9XS  2HIX,11X,2HIY,6X,2HCR,12X,5HSIGMA) 

118  WRITE (6, 119 t THETA, (KK l NM) , NM= 1, 6 ) , CHE K , KX, KY , VR, TENGY 

119  FORMAT  <5X,F5.2,6I2X,F6.3) 2X, F 5 . 3 , 2 (2X , Fll . 4) , 2X , F 8 . 2 , 2X , F U. 4) 

121  CONTINUE 
MX*NNX*i 

C  FOLLOWING  DETERMINES  REPEAT  OF  PRECEEUING  PROCESSES  UNTIL  MAXIMUM 
C  INCIDENT  ANGLE  SPECIFIED  IS  ACHIEVED. 

IF  INNX-NNN)  10  0,1  00,122 

122  CONTINUE 
STOP 
END 
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sea  ice  thickness  using  pulse  echo  techniques  is  demonstrated,  and  recommen¬ 
dations  for  further  study  are  made.  Preliminary  tests  to  measure  Lamb  wave 
propagation  were  successfully  made  and  the  results  are  reported.  (U) 
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